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Abstract 

The  use  of  sodium  laser  beacons  in  adaptive  optics  is  relatively  novel,  and  there  is 
much  current  research  in  this  field.  Perspective  elongation  is  a  side  effect  of  using 
these  beacons,  and  results  in  a  different  shaped  beacon  viewed  by  each  subaperture, 
depending  on  its  location  with  relative  to  the  laser.  Typical  calibration  of  wavefront 
sensors  assumes  that  each  subaperture  image  is  identical,  reducing  accuracy. 

Unfortunately,  telescopes  built  with  a  Coude  path  have  a  rotating  exit  pupil  that 
rotates  when  the  telescope  elevation  or  azimuth  changes.  This  rotation  ensures  that 
static  calibration  of  the  sensor  with  a  single  elongated  reference  beacon  is  inadequate, 
and  a  different  approach  is  required. 

This  research  models  the  sodium  beacon,  its  transmission  through  the  atmosphere 
and  its  measurement  by  a  Shack-Hartmann  wavefront  sensor  (SHWFS).  By  predicting 
the  extent  of  beacon  elongation  and  Coude  rotation,  it  is  possible  to  produce  reference 
images  for  each  subaperture  throughout  an  engagement  scenario.  These  reference 
sources  are  then  used  to  continuously  recalibrate  the  system  as  it  changes  orientation. 

This  model  measures  the  effect  of  perspective  elongation  and  Coude  on  SHWFS 
measurements  to  quantitatively  determine  the  extent  of  degradation  that  occurs.  It 
also  predicts  the  improvement  that  could  be  achieved  by  considering  these  effects, 
and  allows  comparisons  between  alternate  system  configurations.  From  the  computer 
model,  it  was  determined  that  accounting  for  perspective  elongation  and  rotation  can 
reduce  errors  in  measurements  by  up  to  50%  in  a  typical  scenario. 
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EFFECT  OF  COUDE  PUPIL  ROTATION  ON  SODIUM  LASER  BEACON 

PERSPECTIVE  ELONGATION 

I.  Introduction 

Adaptive  Optics  (AO)  is  a  technique  used  to  measure  and  improve  image  quality 
that  has  been  distorted  by  optical  aberrations.  It  is  primarily  used  in  astronomy, 
where  the  atmosphere  is  the  cause  of  aberrations,  and  also  in  retinal  imaging.  A 
natural  star  can  often  be  used  as  a  reference  to  measure  an  atmospheric  effect  on  an 
incoming  wavefront,  but  may  not  be  available  in  the  region  of  interest.  To  compensate 
for  this,  the  sodium  laser  guide  star  was  developed  for  the  Starfire  Optical  Range 
(SOR).  The  laser  beacon  produces  a  light  source  by  resonating  sodium  atoms  in  the 
upper  atmosphere,  and  has  the  advantage  of  being  capable  of  being  moved  over  any 
sky  region  of  interest. 

The  Starfire  Optical  Range  also  operates  a  telescope  that  utilizes  a  Coude  path. 
A  Coude  path  uses  multiple  mirrors  to  reflect  the  incoming  light  in  such  a  way  that 
the  optical  path  beyond  these  mirrors  remains  fixed  in  one  direction,  regardless  of 
telescope  orientation.  Without  using  a  Coude  path,  most  AO  components  are  phys¬ 
ically  attached  to  the  telescope,  severely  limiting  the  allowable  size  and  complexity 
of  optical  components.  A  Coude  path  allows  equipment  to  be  physically  separated 
from  the  telescope,  enabling  greater  flexibility  in  the  design  of  downstream  optics. 
One  unfortunate  effect  of  using  a  Coude  path  is  the  rotation  of  the  optical  image  as 
the  telescope  varies  in  azimuth  or  elevation.  In  many  scenarios,  this  image  rotation  is 
of  little  concern,  and  can  be  corrected  using  post-processing.  If  used  in  conjunction 
with  a  side-launched  laser  beacon,  however,  asymmetries  in  the  reference  source  can 
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degrade  AO  system  performance. 


1.1  Problem  Statement  and  Hypothesis 

The  Starfire  Optical  Range  operates  a  3.5m  telescope  on  a  Coude  path  with  a 
sodium  beacon  reference  for  its  AO  system.  This  combination  produces  the  unique 
effect  of  both  elongating  the  reference  source,  and  also  rotating  its  relative  location 
to  the  sensor.  The  goal  of  this  research  is  to  model  the  behavior  of  Shack- Hartmann 
Wavefront  Sensor  (SHWFS)  performance  with  an  AO  system  using  a  telescope  with 
a  Coude  path  and  a  side-launched  sodium  Laser  Guide  Star  (LGS).  Specifically,  the 
objectives  include: 

•  determine  the  extent  of  Coude  rotation  in  a  typical  engagement  scenario, 

•  model  a  sodium  laser  beacon  and  its  elongation  effects  for  each  subaperture 
within  a  typical  AO  system, 

•  determine  the  effect  of  beacon  elongation  and  Coude  rotation  on  sensor  centroid 
measurements,  and  the  magnitude  of  the  resulting  error  after  reconstructing  of 
the  original  wavefront, 

•  test  out  a  novel  means  of  compensating  for  beacon  elongation  measurement 
errors,  and 

•  maintain  robustness  in  modeling  the  scenario,  to  allow  changes  to  the  engage¬ 
ment  scenario  and  system  configuration  with  minimal  changes  to  source  code. 

1.2  Thesis  Overview 

Chapter  II  provides  the  background  information  required  to  develop  and  under¬ 
stand  the  model  of  this  particular  AO  system  along  with  the  most  relevant  related 
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literature.  Specifically,  it  details  the  cause  and  effect  of  Coude  rotation  and  per¬ 
spective  elongation,  in  addition  to  the  basics  of  AO,  wavefront  measurement  and 
reconstruction.  Chapter  III  details  the  development  of  the  computer  model  and  the 
mathematical  principals  used  to  characterize  the  scenario. 

Chapter  IV  presents  the  results  that  are  produced  by  the  computer  model,  provid¬ 
ing  comparisons  between  the  current  system  configuration  and  the  improvement  that 
could  be  achieved  through  considering  the  effect  of  Coude  rotation  and  perspective 
elongation  during  sensor  calibration. 

Finally,  Ch.  V  draws  conclusions  from  the  results  presented  in  Ch.  IV,  and 
summarizes  the  key  outcomes.  It  also  suggests  how  to  improve  the  current  system, 
and  outlines  areas  that  would  benefit  from  additional  research  in  this  area. 
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II.  Background 


This  chapter  introduces  the  concepts  of  atmospheric  turbulence  and  its  effect  on 
the  propagation  of  light.  Adaptive  optics  (AO)  systems  are  discussed  as  a  method  to 
measure  and  compensate  for  turbulence,  and  sodium  laser  guide  stars  in  particular 
are  discussed.  Finally,  Coude  rotation  is  described  to  enable  a  greater  understanding 
of  the  specific  issues  affecting  the  sponsor  of  this  work.  Current  literature  relating  to 
laser  guide  star  perspective  elongation  is  also  reviewed. 

2.1  Adaptive  Optics 

Light  that  has  been  transmitted  through  a  turbulent  medium  can  be  corrected  by 
an  AO  system  and  returned  to  its  original  form.  An  AO  system  typically  consists  of  a 
Wavefront  Sensor  (WFS)  that  measures  the  turbulence  affecting  the  incident  light,  a 
Fast  Steering  Mirror  (FSM)  that  corrects  low-order  aberrations,  a  Deformable  Mirror 
(DM)  capable  of  correcting  higher-order  aberrations,  and  a  controller  that  commands 
the  mirrors  based  on  WFS  measurements.  The  primary  focus  of  this  research  is  the 
WFS,  and  in  particular  a  proposed  calibration  technique  that  is  designed  to  improve 
the  accuracy  of  WFS  measurements.  The  application  is  for  compensating  atmospheric 
imaging  systems,  so  this  chapter  begins  with  a  review  of  optical  turbulence. 

2.1.1  Atmospheric  Turbulence. 

The  atmosphere  is  not  a  uniform  medium.  Localized  inhomogeneities  in  the  atmo¬ 
sphere  cause  differences  in  temperature,  pressure  and  consequently  index  of  refraction. 
These  localized  inhomogeneities  are  called  eddies,  and  are  defined  as  the  volume  over 
which  the  refractive  index  remains  relatively  constant.  The  size  of  these  eddies  can 
vary,  but  the  minimum  and  maximum  size  of  an  eddy  is  defined  by  the  inner  and 
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Various  PSD  Models 


outer  scales,  Iq  and  L0,  respectively.  Kolmogorov  characterized  the  Power  Spectral 
Density  (PSD)  of  atmospheric  turbulence  <f>  as  a  function  of  wavenumber  k  (measured 
in  [rad/m])  by [1] 

$(«)  =  0.033C^-11/3,  (1) 

where  is  the  structure  constant  of  the  variation  of  the  refractive  index  in  the 
atmosphere  and  effectively  describes  the  strength  of  the  turbulence.  Equation  (1) 
holds  when  1/k  is  within  the  inner  and  outer  scales,  or  alternatively, 

11 

(2) 

-Do  UJ 

Wavenumbers  that  meet  this  condition  are  said  to  be  within  the  inertial  subrange  of 
the  turbulence.  In  Kolmogorov  turbulence,  the  inner  and  outer  scales  are  set  to  0 
and  infinity,  respectively.  These  values  are  not  physically  realistic,  and  more  realistic 
PSD  models  have  been  developed  to  account  for  this,  including  the  Tatarskii,  von 
Karman  and  modified  von  Karman  spectra[20].  These  models  contain  additional 
factors  accounting  for  inner  and  outer  scale  that  are  responsible  for  the  low-frequency 
roll  off  and  the  high  frequency  cut-off,  as  shown  in  Fig.  1. 
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Figure  2.  The  effect  of  eddies  on  an  incoming  wavefront. 

For  light  that  propagates  through  turbulence,  the  different  indices  of  refraction 
within  various  eddies  cause  the  light  to  travel  at  different  speeds.  Figure  2  illustrates 
this  by  showing  various  wavefronts  of  light  as  it  is  transmitted  from  a  star  to  earth. 
Although  a  star  emits  light  spherically,  it  is  approximately  planar  by  the  time  it 
reaches  the  atmosphere  of  earth.  As  it  passes  through  the  atmosphere,  the  plane 
wave  passes  through  eddies  with  varying  refractive  indices,  causing  the  light  to  refract 
differently  at  different  locations.  By  the  time  it  reaches  a  ground-based  telescope,  the 
wavefront  is  no  longer  planar,  but  aberrated. 

Detailed  analysis  of  weak  turbulence  effects  on  wave  propagation  uses  the  Rytov 
method[l].  This  method  essentially  solves  Maxwell’s  equations  with  perturbation 
theory.  The  statistics  of  the  perturbations  are  governed  by  Kolmogorov  theory,  and 
solutions  are  obtained  for  the  turbulence-degraded  field’s  mutual  coherence  function 
and  similar  quantities.  The  details  are  not  presented  here,  just  the  key  resulting 
turbulence  parameters. 

An  aberrated  wavefront  can  be  divided  into  segments  that  are  roughly  spatially 
coherent,  and  this  scale  scale  size  indicates  the  magnitude  of  the  turbulent  phase 
fluctuations.  The  length  of  this  segment  is  known  as  Fried’s  parameter,  tq.  In  the 
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case  of  an  initially  planar  wavefront,  a  more  formal  definition  for  r0  is [9] 


r0  =  0.185 


cum 


(3) 


where  A  is  the  mean  wavelength  of  the  incoming  light  and  JQ~  d £  is  the  integral  of 
the  path  though  which  the  light  is  propagating.  A  small  r0  indicates  a  small  region 
of  spatial  coherence,  and  therefore  large  phase  fluctuations.  Conversely,  a  large  r$ 
indicates  weaker  phase  fluctuations.  For  light  at  500nm  propagating  vertically,  typical 
values  of  r0  range  from  5-10cm. 

In  addition  to  eddies  causing  spatial  variations  in  the  refractive  index  of  the  at¬ 
mosphere,  they  do  not  remain  stationary,  and  hence  cause  temporal  temporal  fluc¬ 
tuations  as  well.  One  parameter  that  is  commonly  used  to  determine  the  rate  of 
temporal  fluctuations  is  the  Greenwood  frequency,  fa-  Typical  values  of  fa  are  of  the 
order  of  10-100Hz[l],  A  formal  definition  is  not  provided  here  because  the  dynamic 
aspect  of  turbulence  is  not  considered  in  this  research. 


2.1.2  Wavefront  Measurement. 

An  aberrated  wavefront  can  be  measured  by  a  wavefront  sensor  (WFS),  which  is 
required  to  compensate  for  the  effects  of  atmospheric  turbulence.  One  common  WFS 
is  the  Shack-Hartmann  (SH)  WFS.  The  SHWFS  comprises  of  an  array  of  lenslets  of 
size  on  the  order  of  r'o .  This  ensures  that  the  light  incident  on  each  lenslet  is  roughly 
spatially  coherent,  and  can  therefore  be  represented  as  a  plane  wave.  As  the  plane 
wave  passes  through  the  lens,  it  is  focused  to  a  sensor  in  the  focal  plane  behind  it. 
Figure  3  depicts  a  lens  focusing  light  onto  a  detector,  and  represents  one  lenslet  within 
a  SHWFS. 

The  light  is  not  be  imaged  to  a  perfect  point,  as  predicted  by  geometric  optics, 
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Figure  3.  The  effect  of  a  lenslet  on  an  incoming  wavefront. 


but  instead  is  affected  by  the  finite  dimensions  of  the  lenslet.  The  resulting  diffraction 
limited  spot  is  the  Fraunhofer  diffraction  pattern  of  the  lenslet,  as  described  by  Ref. 


[10] 


h(u,v)  = 

\Zi 


2n 

P(x,  y)  exp  — j  — — -  (ux  +  vy)  ^  dxdy, 


(4) 


where  A  is  the  amplitude  of  the  incident  field,  A  is  the  wavelength  of  the  incident 
light,  Zi  is  the  propagation  distance  and  P(x,y )  is  the  pupil  function.  For  a  square 
lenslet,  Eq.  (4)  simplifies  to 


h(u'v)  =  s; 

Ad 2 


d/2 


exp 


-d/2 


2n 

-j  — — -  (ux  +  vy)  [>  dxdy, 


A  Zi 


sinc(du)sinc(dv ), 


(5) 

(6) 


where  d  is  the  width  of  the  lenslet. 

Equation  (6)  describes  the  coherent  impulse  response  of  the  system,  which  can  be 
used  to  determine  the  Point  Spread  Function  (PSF)  of  the  system,  Sj(u,  v )  =  | h(u,  v)  |2. 
The  PSF  is  quite  useful  in  incoherent  imaging  systems  because  it  can  be  used  to 
determine  the  image  of  any  input,  not  just  a  point  source.  For  other  inputs,  the 
diffraction  image  Ii(u,  v )  is  simply  a  convolution  of  the  geometrically  predicted  image, 


Figure  4.  Quad  cell  detector. 


Ig(u,v )  and  the  PSF, 


Ii(u,v)  =  Ig(u,v)  ®fi(u,v). 


(7) 


The  light  focussed  by  a  Shack- Hartmann  lenslet  is  measured  by  an  array  of  pixels 
at  the  focal  plane  of  the  lenslet.  Figure  4  depicts  the  simplest  detector,  a  quad  cell 
detector  consisting  of  a  2x2  pixel  array.  The  individual  detector  pixels  are  referenced 
according  to  their  vertical  and  horizontal  location  relative  to  the  center  of  the  detector. 
Because  the  imaged  spot  typically  spans  multiple  pixels,  the  location  of  its  center  can 
be  determined  from  the  relative  number  of  photons  measured  at  each  pixel.  The  x 
and  y  centroids,  xc  and  yc,  can  be  computed  by  the  weighted  centroid  calculation 


xc 


Vc 


5 

i,j  lij 


(8) 


where  Xij  and  ijij  are  the  numeric  indices  of  each  pixel,  indicating  their  location  with 
respect  to  the  center  of  the  detector,  and  is  the  measured  intensity  in  the  i,j 
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pixel.  For  the  quad  cell  detector  depicted  in  Fig.  4,  this  simplifies  to 


Ib  +  Ip  —  (I A  +  Ic) 

Ia  +  Ib  +  Ic  +  Id 

IA  +  Ib  ~  ( Ic  +  Ip)  /q\ 

Ia  +  Ib  +  Ic  +  Id  ’  U 

where  Ia,b,c,d  represents  the  measured  irradiance  in  pixels  A,B,C  and  D. 

For  a  plane  wave  perpendicular  to  the  optical  axis,  the  image  is  formed  on  the 
optical  axis.  For  plane  waves  oblique  to  the  optical  axis,  however,  the  image  is 
focused  to  a  point  that  is  laterally  displaced  from  the  optical  axis,  as  illustrated  in 
Fig.  3.  Although  Fig.  3  shows  only  a  one  dimensional  cross  section  of  the  light  being 
imaged  by  the  lens,  the  light  is  actually  displaced  in  both  the  horizontal  and  vertical 
directions,  depending  on  the  x  and  y  slopes  of  the  incident  plane  wave.  In  a  SHWFS, 
a  detector  is  placed  at  the  focal  length  of  the  lenslet  and  measures  the  intensity  of 
the  light  in  each  pixel  of  the  detector. 

From  these  centroid  measurements,  the  original  x  and  y  slopes  of  the  incident 
wavefronts  can  be  determined  using  the  geometry  featured  in  Fig.  3.  Assuming  the 
paraxial  (small  angle)  approximation  such  that 

sin(0)  «  9,  (10) 

the  slope  of  the  incident  wavefront,  9w ,  can  be  found  from  simply  the  lens  focal 
length,  fi  and  the  spot  displacement,  ds  shown  in  the  Fig.  3.  Using  trigonometry, 

sin(6»T)  =  ^  «  9t,  (11) 

Ji 

and  through  similar  triangles, 

9W  =  0t,  (12) 


xc  = 

Vc  = 
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Figure  5.  Deformable  mirror  correcting  the  aberrations  of  an  incident  wavefront. 


such  that 


a  ^ s 

Ji 


(13) 


Since  the  displacement  of  the  spot  is  therefore  determined  by  the  slope  of  the 
incident  wavefront,  the  centroid  theoretically  measure  the  intensity- weighted  phase 
gradient.  The  measured  x  and  y  centroids  are  used  to  find  the  x  and  y  slopes,  Sx  and 
Sy.  From  Eq.  (13), 


c  %c 

Sx  ~  i 
Ji 

Sy  ~  (14) 

2.1.3  Wavefront  Reconstruction. 

The  measured  slopes  from  a  wavefront  sensor  provide  enough  information  to  re¬ 
construct  the  shape  of  the  incident  wavefront.  Since  the  reconstructed  wavefront  is 
the  phase  of  a  wavefront  that  was  originally  planar,  it  also  represents  the  Optical 
Path  Difference  (OPD)  of  the  wavefront  for  different  points  across  the  entrance  pupil. 
A  segmented  or  deformable  mirror  can  be  applied  to  the  light  to  compensate  the 
OPD  across  the  wavefront.  The  resulting  ‘corrected’  image  is  much  less  aberrated,  as 
illustrated  in  Fig.  5. 
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4Vm+l)  ^(n+l  ,m+l) 
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Figure  6.  The  Fried  geometry,  as  defined  by  the  location  of  actuators  with  respect  to 
sensor  measurements. 

The  location  at  which  DM  actuators  are  placed  with  respect  to  the  SHWFS  sub¬ 
apertures  is  called  the  system  geometry.  The  ‘Fried  geometry’  is  defined  by  actuators 
placed  at  the  corners  of  a  subaperture,  rather  than  directly  over  the  center,  where 
the  slope  measurements  are  made.  This  is  illustrated  in  Fig.  6.  In  this  figure,  the 
n,m  subaperture  is  shown,  with  x  and  y  slopes  Sx  and  Sy,  and  reconstructed  phase 
measurements  (f)n,m  that  are  to  be  applied  to  the  DM. 

The  relationship  between  slope  measurements  and  reconstructed  phase  can  be 
described  using  linear  algebra,  where 

Sx(n,m)  g  (077+1,777  077, 777)  T  g  (071+1,771+1  0?7,777+l) 

Sy(n,m)  g  (07i,?7i+l  077, 771)  T  (071+1,771+1  077+1,777)  (15) 

These  equations  can  be  written  in  matrix-vector  form,  such  that 

s  =  G  0.  (16) 
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Unfortunately,  G  is  not  a  square  matrix,  so  the  equation  cannot  be  solved  for  <fi  using 
simple  matrix  inverse  techniques.  Instead,  a  pseudo-inverse  method  can  be  used  to 
retrieve  phase  values  according  to  [11]. 


(f>  =  G+s 

4>  =  (GtG)-1  Grs 

4>  =  Hs,  (17) 


where  H  is  the  reconstructor  matrix. 

Typically,  approximately  87%  of  the  phase  variance  caused  by  atmospheric  turbu¬ 
lence  can  be  represented  by  tilt  [11],  and  can  be  corrected  using  a  rigid  mirror.  This 
tilt  is  removed  by  a  FSM,  which  greatly  reduces  the  requirements  placed  on  the  DM. 
Additionally,  it  can  typically  be  updated  much  faster  than  the  DM,  using  basic  tilt 
measurements  [such  as  Eq.  (9)]  directly  from  a  tracking  camera,  rather  than  using 
data  from  the  WFS.  A  diagram  of  a  standard  AO  system  using  all  these  components 
is  shown  in  Fig.  7. 

2.2  Laser  beacons 

The  SHWFS  described  above  requires  a  reference  source  to  measure  the  effects  of 
the  atmosphere  on  a  known  wavefront.  Traditionally,  stars  have  been  used  as  reference 
sources  because  their  distance  provides  an  almost  planar  wave  by  the  time  the  light 
reaches  Earth,  despite  originally  being  emitted  spherically.  This  is  illustrated  in  Fig. 
8. 

Unfortunately,  natural  guide  stars  suffer  from  two  major  drawbacks.  Firstly,  stars 
produce  a  relatively  small  amount  of  light  to  be  used  as  a  reference  source,  making  the 
measurement  more  susceptible  to  noise.  More  importantly,  a  viable  star  may  not  be 
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Sodium  Layer 


Figure  8. 
propagate. 


Spherical  light  waves  emitted  from  a  star  appearing  more  planar  as  they 
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available  in  the  region  of  interest,  especially  when  moving  targets  are  considered.  For 
these  reasons,  artificial  (laser)  guide  stars  are  becoming  more  common  for  adaptive 
optics  uses. 

There  are  two  common  types  of  LGS’s,  Rayleigh  and  sodium.  Rayleigh  beacons 
rely  on  Rayleigh  scattering  from  air  molecules  at  low  altitudes  (approx.  20km)  [19]. 
These  beacons  can  be  produced  for  lasers  of  varying  wavelengths,  allowing  relatively 
cheap  production  of  high  power  guide  stars.  Unfortunately,  while  the  majority  of 
atmospheric  turbulence  occurs  at  lower  altitudes,  the  cone  of  light  measured  from 
Rayleigh  beacons  does  not  include  the  majority  of  the  atmosphere.  Alternatively, 
sodium  LGS’s  use  resonant  scattering  from  sodium  atoms  at  mesospheric  altitudes, 
much  higher  in  the  atmosphere.  The  light  from  sodium  LGS’s  is  therefore  affected 
by  a  much  greater  volume  of  the  atmosphere,  and  they  provide  a  much  more  realistic 
representation  of  the  effect  of  the  total  atmosphere.  This  is  known  as  the  cone  effect, 
and  is  illustrated  in  Fig.  9. 

One  disadvantage  of  sodium  LGS  use  is  the  requirement  for  a  specific  wavelength 
laser.  Figure  10  shows  both  Rayleigh  and  Sodium  layer  scattering  from  a  sodium 
laser.  The  much  longer  streak  is  the  Rayleigh  beacon,  caused  by  low-altitude  Rayleigh 
scattering.  The  much  smaller  spot  in  the  upper  right  of  the  image  (and  inset)  is  the 
sodium  guide  star. 

2.2.1  Sodium  Layer. 

In  the  upper  atmosphere  at  an  altitude  of  approximately  90km,  there  is  a  layer  of 
sodium  atoms  at  a  relatively  high  concentration [19].  While  the  origin  of  these  atoms 
is  not  certain,  research[16]  suggests  that  it  may  be  caused  by  meteorites  boiling  off 
sodium  as  they  pass  through  the  atmosphere.  The  exact  height  and  depth  of  this  layer 
varies  with  many  factors,  including  time  and  location,  but  the  literature  commonly 
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Rayleigh  LGS 


Figure  9.  The  cone  effect.  Due  to  the  higher  altitude  of  the  Sodium  LGS,  a  greater 
volume  of  the  atmosphere  affects  the  LGS,  providing  a  more  accurate  reference  beacon. 


Figure  10.  Sodium  laser  guide  star  showing  Rayleigh  scattering  and  sodium  layer 
Resonance.  Photo  by  Art  Goodman  and  inset  by  Jack  Drummond. 
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places  the  layer  between  the  altitudes  of  approximately  85km-95km[3,  19]  at  zenith. 
The  concentration  of  sodium  atoms  within  the  layer  typically  follows  a  Gaussian 
distribution  as  a  function  of  altitude[16]. 

2.2.2  Sodium  Laser. 

A  laser  operating  at  the  589. 2nm  sodium  wavelength  can  be  created  by  using  a 
Sum-Frequency  Generation  (SFG)  technique.  In  this  technique,  two  NckYAG  lasers, 
one  operating  at  1064  nm  and  another  at  1319  nm  can  be  combined  to  produce  a 
laser  of  wavelength  589  nm[14],  SOR  has  previously  demonstrated  laser  output  levels 
at  the  sodium  wavelength  of  approximately  50  W[8]. 

2.2.3  Laser  Guide  Star  Perspective  Elongation. 

One  of  the  issues  affecting  the  performance  of  LGS  systems  is  perspective  elonga¬ 
tion.  This  is  caused  by  the  finite  depth  of  the  sodium  layer  and  laser  beam,  which 
interact  to  form  a  ‘column’  of  light  in  the  upper  atmosphere.  When  viewed  directly 
below  the  column,  the  scattering  spot  appears  circular,  similar  to  how  a  star  might 
appear.  When  viewed  off-axis,  however,  the  full  height  of  the  column  becomes  visi¬ 
ble,  and  the  spot  becomes  less  circular  and  more  elliptical  in  appearance.  The  further 
off-axis,  the  greater  the  elongation  of  the  imaged  spot.  This  effect  can  be  witnessed 
in  Fig.  11. 

This  issue  is  particularly  important  to  consider  when  using  a  SHWFS.  Centroid 
measurements  become  less  accurate  when  modeled  with  an  ideal  spot  and  not  the 
more  realistic  elliptical  shape.  The  geometric  details  are  shown  in  Figs.  13-15  and 
described  in  Sec.  3.1.3.  These  measurements  affect  the  wavefront  slope  calculations, 
and  consequently  the  reconstructed  phase  of  the  wavefront  [21], 
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Figure  11.  Laser  beacon  elongation  from  measurements  taken  at  the  Keck  observatory 
in  Hawaii.  In  this  image,  the  laser  is  launched  from  the  right  side  of  the  sensor.  [26]. 

2.3  Coude  Rotation 

Unlike  many  telescopes,  telescopes  that  use  a  Coude  path  are  designed  to  provide 
a  fixed  optical  axis  regardless  of  telescope  elevation  and  azimuth.  Unfortunately,  this 
causes  the  optical  image  to  rotate  as  the  telescope  changes  its  orientation.  A  de¬ 
rotation  system  can  be  used  to  compensate  for  this  effect  and  therefore  maintain  a 
pupil  and  image  orientation  that  remains  fixed,  but  the  additional  components  reduce 
the  optical  throughput  to  the  wavefront  sensor  [17].  Omitting  the  de-rotator  creates 
additional  complications  when  using  laser  beacons,  due  to  the  effects  of  laser  beacon 
elongation,  as  each  subaperture  spot  shape  is  different,  as  shown  in  Fig.  11.  When 
the  telescope  moves,  the  pupil  and  laser  beacon  launch  position  rotate  relative  to  the 
wavefront  sensor,  so  these  dynamically  change  shape. 

From  Ref.  [13],  the  rotation  angle  from  a  Coude  telescope  is 

@Coude  @ Gregory  A  @ el  @az  A  C  (18) 
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where  9ei  is  the  telescope  elevation  angle,  9az  is  the  telescope  azimuth  angle,  C  is 
a  constant  that  describes  the  initial  relative  rotation  of  sensors,  and  9 Gregory  is  the 
orientation  of  the  object  at  Gregorian  focus.  For  a  stationary  object  such  as  a  laser 
beacon,  there  is  no  rotational  component  at  the  Gregorian  focus,  such  that  9oregory  = 
0.  From  the  perspective  of  a  fixed  detector,  the  relative  angle  of  the  laser  beacon  with 
respect  to  the  sensor  is  simply  the  negative  of  the  Coude  rotation  of  the  exit  pupil 
with  respect  to  a  fixed  telescope,  so  that 


f)  fj 

ubeac  u  Coude 

@beac  @az  @el  C . 


(19) 


Equation  (19)  shows  that  altering  the  elevation  or  azimuth  of  the  telescope  rotates 
the  relative  location  of  the  laser  beacon,  which  correspondingly  affects  the  elongation 
of  LGS  spots  on  individual  detectors  as  the  telescope  moves. 

2.4  Starfire  Optical  Range 

The  Starfire  Optical  Range  is  located  at  the  Kirtland  Air  Force  Base  (AFB)  in 
Albuquerque,  New  Mexico.  It  is  a  laboratory  facility  that  forms  part  of  the  Directed 
Energy  Directorate  at  the  Air  Force  Research  Laboratory  (AFRL).  SOR  houses  a 
handful  of  telescopes  capable  of  tracking  earth-orbiting  satellites  [28].  The  most  no¬ 
table  telescope  is  3.5m  in  diameter  and  is  shown  with  its  sodium  laser  beacon  in  Fig. 
12.  The  telescope  forms  part  of  an  adaptive  optics  system  capable  of  compensating 
for  the  effects  of  atmospheric  turbulence. 

SOR  also  utilizes  a  side-launched  sodium  laser  guide  star,  and  was  the  first  site 
to  demonstrate  closed-loop  laser  beacon  adaptive  optics[ll].  Updates  planned  for 
the  range  in  2009  propose  a  SHWFS  using  a  32x32  array  of  quad  cell  detectors [15]. 
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Figure  12.  3.5m  telescope  at  the  Starfire  Optical  Range  near  Albuquerque,  NM.  with 
side  lauched  laser[27]. 

It  was  predicted  that  the  unvignetted  Field  of  View  (FOV)  of  the  system  would  be 
approximately  0.5  mrad  (100  arcsec). 

2.5  Related  Research 

This  section  outlines  the  recent  research  relating  to  the  compensation  of  laser 
beacon  elongation  effects  on  various  systems.  Several  proposals  that  directly  relate 
to  improving  measurements  affected  by  perspective  elongation  are  described,  usually 
requiring  expensive  or  complex  solutions.  The  research  described  in  Chapters  III-V 
investigates  the  feasibility  of  implementing  a  much  lower  risk  solution  than  can  be 
implemented  in  software  alone. 

2.5.1  Improved  Noise  Model  with  Multiple  Sensors. 

Reference  [2]  reports  the  findings  of  a  closed-loop  ground-layer  AO  system  using 
laser  guide  stars,  specifically  for  an  Extremely  Large  Telescope  (ELT).  Since  the  extent 
of  beacon  elongation  is  dependent  on  the  ground  separation  between  each  subaperture 
and  the  laser  launch  position,  it  is  a  major  concern  for  very  large  telescopes.  Even  in 
the  preferred  center-launched  laser  configuration,  their  is  still  a  separation  between 
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laser  and  subaperture  of  up  to  21m  for  a  42m  diameter  ELT,  compared  with  the 
maximum  separation  of  4.5m  currently  possible  at  SOR. 

Reference  [2]  outlines  the  development  of  an  improved  noise  model  to  account  for 
beacon  elongation  in  their  simulations.  It  describes  how  the  variance  in  measurement 
error  is  greatest  in  the  axis  of  elongation,  and  how  using  multiple  side-launched  lasers 
with  corresponding  WFS’s  can  perform  almost  as  well  as  a  center  launched  laser,  if 
this  noise  model  is  considered.  They  concede  that  the  validity  of  this  model  requires 
further  study. 

2.5.2  Pulsed  Laser. 

Reference  [4]  offers  a  unique  solution  to  compensate  for  perspective  elongation 
caused  by  the  thickness  of  the  sodium  layer.  By  substituting  a  narrow  pulsed  laser 
in  lieu  of  a  Continuous  Wave  (CW)  laser,  they  propose  that  the  depth  of  the  sodium 
layer  that  is  illuminated  by  the  laser  at  any  time  can  be  shallow  enough  that  minimal 
elongation  of  the  LGS  is  incurred. 

Unfortunately,  this  solution  has  complications.  Although  the  proposed  pulsed 
lasers  exist  (requiring  pulse  lengths  <2ps  and  enough  power  for  adequate  photon 
return),  the  WFS  itself  must  be  upgraded  to  work  effectively  with  the  laser.  It  is 
proposed  that  either  the  image  is  realigned  or  the  WFS  is  moved  as  the  illuminated 
volume  progresses  through  the  sodium  layer.  While  technically  possibly,  this  solution 
still  requires  significant  work  in  sensor  design  and  implementation. 

2.5.3  Matched  Filtering. 

Reference  [22]  describes  several  techniques  designed  to  improve  SHWFS  perfor¬ 
mance,  including  threshold  center  of  gravity  and  weighted  center  of  gravity.  Although 
the  techniques  were  proposed  to  analyze  noise  susceptibility  using  a  symmetric  LGS, 
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one  technique  may  prove  to  be  well-suited  to  compensating  perspective  elongation 
effects.  Reference  [22]  proposes  image  correlation  as  a  form  of  matched  filtering  to 
determine  the  actual  centroid  of  a  detected  spot.  The  method  involves  correlating 
an  ideal  image  with  the  actual  detected  image  to  determine  the  location  of  greatest 
similarity. 

This  method  may  severely  reduce  measurement  errors  caused  by  elongation,  but 
suffers  two  major  drawbacks.  First,  it  would  require  ideal  images  for  each  subaperture 
and  each  telescope  orientation,  which  requires  more  research  into  LGS  production  or 
experimental  data.  Secondly,  the  correlation  operation  is  computationally  intensive, 
and  may  not  be  a  feasible  solution  for  real-time  operation  with  today’s  processing 
technology. 

2.5.4  Centroid  Gain  Correction. 

Reference  [25]  describes  simulations  designed  at  SOR  that  model  the  effect  of 
perspective  elongation  on  centroid  gain  measurements.  The  research  is  designed  to 
consider  a  dynamic  laser  launch  position  (caused  by  Coude  rotation).  The  simulations 
include  a  Gaussian  beacon  model  for  subaperture  images  which  is  used  to  determine 
centroid  gains  at  each  subaperture  by  performing  a  calibration  type  operation  with 
each  image.  The  model  shows  the  greatest  variance  in  gradient  measurements  cor¬ 
responds  to  the  subapertures  affected  by  the  greatest  elongation.  The  model  also 
predicts  the  Strehl  ratio  for  various  scenarios,  predicting  increased  system  perfor¬ 
mance  at  higher  elevation  angles  and  in  weaker  turbulence. 

The  research  described  in  Chapters  III-V  extends  upon  what  is  developed  in  Ref. 
[25]  to  consider  the  effect  of  elongation  after  wavefront  reconstruction.  The  method 
employed  to  determine  centroid  gains  in  Ref.  [25]  is  used  to  predict  new  calibration 
factors  to  improve  system  performance.  The  research  in  this  thesis  goes  further  to 
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investigate  the  effect  of  detector  resolution  on  system  performance. 
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III.  Methodology 


This  chapter  outlines  the  methods  and  techniques  used  to  achieve  the  objectives 
outlined  in  Chapter  1.  It  outlines  the  parameters  used  to  represent  the  SOR  site  and 
the  mathematical  expressions  that  are  used  to  model  the  system  geometry.  Finally, 
the  implementation  of  these  functions  into  MATLAB  environment  is  described. 

3.1  Modeling  the  Site 

As  outlined  in  Chapter  1,  an  objective  of  this  model  is  to  maintain  robustness 
such  that  many  different  AO  system  configurations  may  be  examined  with  minimal 
alteration  of  source  code.  Of  primary  concern,  however,  is  the  current  configuration  in 
operation  at  the  Starfire  Optical  Range  (SOR).  For  this  reason,  the  SOR  configuration 
is  considered  the  baseline  against  which  other  setups  are  compared. 

3.1.1  SOR  Configuration. 

Reference  [15]  discusses  a  proposed  upgrade  to  SOR’s  AO  system  for  its  3.5m 
telescope  in  2009.  The  proposed  upgrade  is  used  as  the  baseline  for  modeling  the  AO 
system,  and  is  described  below: 

•  The  primary  mirror  is  3.5m  in  diameter; 

•  The  WFS  contains  an  array  of  32x32  subapertures;  and 

•  Each  subaperture  is  imaged  on  to  a  quad  cell  detector. 

The  following  assumption  is  made  about  specifications  not  found  in  the  literature: 

•  The  sodium  laser  is  launched  approximately  1  m  from  the  edge  of  the  telescope 
aperture,  based  on  pictures  of  the  telescope  and  laser  launch. 
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3.1.2  Turbulence  Model. 


The  turbulence  model  used  in  this  research  is  Kolmogorov  turbulence,  as  discussed 
in  Chapter  2.  To  simulate  turbulence  for  a  telescope  viewing  directly  overhead,  a  sin¬ 
gle  phase  screen  may  be  produced  that  mimics  the  phase  difference  from  an  originally 
planar  wave  that  atmospheric  turbulence  would  create [20].  To  assist  in  simulating 
wave  propagation,  two  MATLAB  toolboxes  were  used.  These  are  WaveProp  and 
AOTools  toolboxes,  developed  by  the  Optical  Sciences  Company  (tOSC)  [5,  6].  The 
included  function  kolmogphzscrn  can  produce  a  randomly  drawn  phase  screen  to 
model  turbulence  for  a  given  r0  and  screen  size. 

To  determine  the  required  r0  for  the  SOR  site,  it  is  necessary  to  determine  the 

profile  and  wind  model  along  the  propagation  path.  The  Hufnagel- Valley  5/7 
(H-V5/7)  structure  constant  profile  is  commonly  used,  as  is  the  Greenwood/Gaussian 
wind  model,  which  depends  on  parameters  such  as  wind  velocity  and  height  of  the 
tropopause.  A  specific  case  of  the  Greenwood/Gaussian  model  is  the  Button  model [1], 
which  has  been  implemented  in  these  calculations.  Researchers  at  tOSC  have  indi¬ 
cated  that  these  models  give  good  agreement  with  atmospheric  measurements  made 
at  SOR, 

3.1.3  Sodium  Layer  Dimensions. 

The  distance  to  the  sodium  layer  boundaries  is  a  consideration  that  affects  the 
extent  of  beacon  elongation.  Common  values  for  the  height  of  the  sodium  layer  at 
zenith  are  around  85-95km  [24,  25].  Regardless  of  the  height  chosen,  the  distance  to 
the  sodium  layer  increases  as  the  telescope  moves  away  from  the  zenith  position.  The 
geometry  used  to  determine  the  distance  to  the  sodium  layer  boundaries  is  shown  in 
Fig  13. 

The  greater  distance  to  the  sodium  layer  boundaries  at  lower  elevations  typically 
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Figure  13.  Geometry  showing  the  dependency  of  telescope  altitude  on  the  distance  to 
the  upper  boundary  of  the  sodium  layer.  The  blue  lines  represent  the  upper  and  lower 
boundaries  of  the  sodium  layer  surrounding  a  section  of  the  earth  in  green.  6ei  is  the 
elevation  angle  of  the  telescope  above  the  horizon. 

produce  a  smaller  angular  width  of  the  laser  beacon,  corresponding  to  less  elongation. 
This  is  despite  the  longer  column  of  light  produced  by  the  laser  as  it  travels  at  an 
angle  through  the  layer.  Using  the  law  of  cosines,  the  distance  to  the  (upper  or  lower 
edge)  sodium  layer  boundary  for  a  given  elevation  angle,  9ei,  is 

hNa  —  R  cos  (^2  +  +  \j (hzen  +  R)  —  R2  sin2  ^—  +  9ei ^ . ,  (20) 

where  R  =  6371  km  is  the  average  radius  of  the  earth  and  hzen  is  the  (upper  or  lower) 
height  of  the  sodium  layer  at  zenith.  At  zero  elevation,  the  distances  to  the  lower  and 
upper  boundaries  of  the  sodium  layer  are  therefore  104.4  and  110.4  km,  respectively. 

3.1.4  Angular  Width. 

For  side  mounted  lasers,  such  as  what  is  operated  at  SOR,  the  laser  must  be 
tilted  to  center  the  beacon  over  the  telescope.  This  amount  of  tilt  required  affects  the 


26 


^upper 


Figure  14.  Geometry  showing  the  effect  of  laser  tilt  and  distance  to  the  sodium  layer 
boundaries  on  the  angular  width  of  the  imaged  beacon  at  a  particular  subaperture [11]. 


angular  width  of  the  beacon  as  viewed  from  a  ground  position.  Figure  14  illustrates 
how  the  laser  tilt  and  distance  to  the  sodium  layer  boundaries  are  related  to  the 
angular  width  of  the  imaged  beacon. 

Applying  trigonometry  to  Fig.  14,  the  angular  width  of  the  beacon,  relative  to  a 
subaperture  within  the  telescope  is 
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where  Dupper  and  Diower  are  projected  horizontal  distances  from  the  subaperture  to 
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Figure  15.  Illustration  of  the  telescope  aperture  showing  relative  distance  from  each 
subaperture  to  the  laser  beacon. 

the  edge  of  the  scattering  spot,  as  shown  in  Fig.  14.  Dupper  and  Diower  are  given  by 

DUpper  Drei  h]ya_u  sill  (T/?//)  , 

DUpper  Drei  h]ya_u  sill  (dtut)  , 

Drei  is  the  distance  from  each  subaperture  to  the  laser,  and  h^a _u  and  h^a_i  are  the 
distances  to  the  upper  and  lower  boundaries  of  the  sodium  layer,  respectively.  Drei 
depends  on  the  rotation  of  the  laser  beacon  relative  to  the  telescope,  as  calculated 
from  Eq.  (19),  and  the  location  of  each  subaperture  within  the  telescope.  This  is 
illustrated  in  Fig.  15. 

3.1.5  Satellite  Track. 

To  compare  the  performance  of  various  system  configurations,  two  sample  satellite 
paths  have  been  analyzed.  The  first  is  the  International  Space  Station  (ISS),  as  it 
could  have  been  viewed  while  it  passed  Albuquerque,  NM  on  Wednesday,  October 


Figure  16.  Track  of  two  satellites  as  they  pass  over  Albuquerque,  NM  on  20  October, 
2010.  The  blue  track  is  the  LEO  International  Space  Station,  orbiting  from  right  to 
left  on  these  axes.  The  ISS  passes  in  under  10  minutes.  The  red  track  is  the  MEO 
GPS  satellite  G28,  orbiting  from  left  to  right  and  taking  7  hours  to  complete  the  pass. 

20,  2010  at  approximately  2130h  local  time.  The  track  data  is  illustrated  in  Fig. 
16,  and  was  obtained  from  the  NASA  website [18].  The  ISS  maintains  a  Low  Earth 
Orbit  (LEO),  and  this  particular  viewing  tracks  almost  180°  in  azimuth,  and  reaches 
a  maximum  elevation  of  54.8°  above  the  horizon.  This  entire  pass  lasts  less  than  10 
minutes. 

The  second  is  the  path  of  a  GPS  satellite,  G28,  viewed  from  the  same  location 
between  0900h  and  1600h  on  the  same  day.  This  satellite  is  orbiting  in  a  Medium- 
Earth  Orbit  (MEO),  and  therefore  takes  much  longer  to  complete  its  trajectory.  This 
viewing  tracks  258°  in  azimuth,  and  reaches  a  maximum  elevation  of  64°  above  the 
horizon.  This  satellite  takes  7  hours  to  complete  this  pass.  The  track  data  was 
extracted  from  satellite  visibility  software  freely  available  on  the  internet  [23]. 

These  satellites  were  chosen  primarily  due  to  the  availability  of  track  data  and 
difference  of  orbit  types.  The  key  benefits  of  these  orbits  are  their  relatively  high 
elevation  angle  and  wide  azimuthal  change  compared  to  other  crossings.  The  greater 
change  in  telescope  orientation  required  to  view  these  satellites  ensures  that  the  Coude 
effect  is  more  evident  than  in  tracks  which  vary  less  in  azimuth  and  elevation.  Eight 


29 


equally  spaced  points  along  each  track  are  analyzed  in  this  research. 

3.2  Modeling  the  Wavefront  Sensor 

The  preceding  section  described  the  environment  being  modeled,  including  the 
physical  characteristics  of  the  target,  transmission  path  and  receiving  telescope.  This 
section  provides  detail  about  how  this  information  is  applied  to  model  the  impact  to 
the  sensor. 

As  in  a  real  AO  system,  the  Shack-Hartmann  sensor  must  be  calibrated  against  a 
known  source  to  determine  the  relationship  between  the  tilt  present  on  a  wavefront 
to  the  resulting  lateral  displacement  of  the  spot  imaged  behind  each  lenslct.  This  is 
typically  achieved  in  a  laboratory  by  using  a  calibration  laser  that  acts  as  a  point 
source,  i.e.,  has  no  finite  extent  nor  elongation. 

Once  wavefront  slopes  have  been  measured  by  the  WFS,  they  are  corrected  ac¬ 
cording  to  this  system  calibration.  The  incident  wavefront  can  then  be  corrected 
through  commands  sent  to  the  DM,  based  on  the  reconstruction  of  these  slope  mea¬ 
surements.  The  research  described  here  uses  this  technique  as  a  baseline  but  also 
models  more  realistic  reference  sources  based  on  the  previous  geometric  discussion. 
It  is  hypothesized  that  these  new  reference  sources  produce  more  accurate  calibration 
factors  and  therefore  reduce  errors  in  the  reconstructed  wavefront. 

Finally,  the  reconstructed  wavefront  is  used  to  quantitatively  determine  how  per¬ 
spective  elongation  and  Coude  rotation  affect  the  performance  of  wavefront  measure¬ 
ments  and  allow  effective  comparisons  between  different  configurations  to  be  per¬ 
formed. 
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3.2.1  Subaperture  Images. 


Modeling  the  unaberrated  shape  that  the  beacon  produces  as  viewed  from  the 
telescope  is  considerably  complex,  due  to  these  factors: 

•  the  laser  beam  is  slightly  diverging  as  it  passes  through  the  sodium  layer, 

•  the  laser  beam  varies  in  spatial  intensity, 

•  the  density  of  the  sodium  atoms  within  the  sodium  layer  varies,  and  depends 
on  many  factors,  and 

•  each  lenslet  views  the  beacon  from  a  unique  perspective  angle,  based  on  its 
location  within  the  aperture. 

Due  to  these  factors,  an  accurate  analytic  model  is  very  difficult  to  achieve.  From 
a  simply  qualitative  perspective,  empirical  data  shows  that  the  beacon  appears  to 
have  an  elliptical  shape,  such  as  in  Fig.  11.  This  model  therefore  uses  a  simple 
two-dimensional  asymmetric  Gaussian  intensity  distribution  to  as  a  first  order  ap¬ 
proximation  of  the  beacon  in  each  subaperture.  This  Gaussian  function  is  given  by 


f(x,y)  =  a  exp 


( X  X0ffset)  ( y  y  off  set) 


2cl 


2cl 


(22) 


where  a  is  a  constant,  x,y  represent  the  spatial  position  in  the  subaperture  with 
respect  to  the  major  and  minor  axes  of  the  ellipse  (in  meters),  x,y0ffset  represent  the 
lateral  displacement  of  the  spot  with  respect  to  the  center  of  the  subaperture,  and 
cx,  cy  represent  the  width  of  the  spot’s  intensity  along  the  major  and  minor  axes  of 
the  ellipse. 

For  simpler  computations,  the  major  and  minor  axes  are  aligned  to  a  Cartesian 
grid,  and  coordinate  transformation  is  used  to  rotate  the  shape  of  each  subaperture 
image  so  that  the  major  axis  is  aligned  with  9iaser.  The  transformation  used  to  create 
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the  rotated  grid  coordinate  system  (x2,  y2 )  from  the  original  coordinates  (xi,  yi)  at 
an  angle  of  9iaser  is  given  by 


X2 

X\  cos  0iaser  y\  sin  0iaser , 

V2 

X\  sin  0iaser  +  y\  COS  0 laser  • 

(23) 

This  transformation  allows  the  angular  width  of  the  minor  axis  to  be  held  constant, 
representing  the  width  of  a  non-elongated  beacon.  For  this  model,  cy  is  kept  at 
0.1,  representing  an  ellipse  about  lm  wide  at  zenith,  which  matches  measurements 
presented  in  Ref.  [25],  and  appears  similar  to  actual  data  presented  in  Ref.  [4], 

As  previously  discussed,  the  angular  width  of  the  major  axis  varies  according  to 
the  perspective  of  the  laser  beacon  as  viewed  from  the  telescope.  To  produce  the 
diffraction  image  present  at  the  detector,  the  Gaussian  function  is  convolved  with 
the  PSF  of  the  Shack-Hartmann  lenslet  in  the  spatial  domain.  The  relative  widths 
of  the  PSF  and  the  Gaussian  functions  therefore  determine  how  much  of  an  effect 
the  elongation  has  on  the  resulting  image.  If  the  width  of  the  beacon  image  is  much 
narrower  than  the  PSF,  the  elongation  is  somewhat  masked  in  the  convolution  process. 
In  the  system  that  is  being  modeled,  the  beacon  image  ranges  from  roughly  the  same 
size  of  the  PSF  to  about  50%  larger,  depending  on  the  location  of  the  subaperture 
and  elevation  of  the  telescope. 

The  convolution  of  the  two  functions  is  performed  using  a  high-resolution  grid  to 
adequately  capture  the  effect  of  the  beacon  shape,  and  the  intensity  of  the  resulting 
diffraction  image  is  integrated  over  each  pixel’s  area  in  the  detector  to  produce  pho¬ 
tocounts.  In  this  model,  the  effect  of  shot  noise  on  the  detected  photocounts  is  not 
considered. 
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Calibration  curves 


Figure  17.  Calibration  curves  for  a  quad-cell  detector.  The  x  axis  shows  the  amount  of 
tilt  (in  radians  of  wavefront  slope)  applied  to  the  imaged  beacon,  and  the  y  axis  shows 
the  location  of  the  measured  centroid  (in  units  of  detector  width).  The  blue  curve 
is  achieved  using  the  basic  sine2  function,  whereas  the  dotted  black  curve  is  from  an 
elongated  beacon  on  the  opposite  side  of  the  telescope  from  the  laser. 


3.2.2  Calibration. 

The  centroid  calculation  in  Eq.  (9)  does  not  perfectly  determine  the  center  of  an 
imaged  spot,  primarily  due  to  the  discrete  nature  of  the  pixels  that  are  measuring 
a  continuous  spot.  The  accuracy  can  be  improved  by  using  more  pixels  within  a 
subaperture,  but  a  basic  quad  cell  is  commonly  used,  as  is  the  case  in  this  model. 

To  account  for  this,  a  calibration  curve  is  often  produced  to  determine  the  mea¬ 
sured  centroid  for  a  given  tilt.  A  source  providing  a  diffraction-limited  spot  is  applied 
with  varying  tilt.  As  the  applied  tilt  is  varied,  the  measured  tilt  typically  follows  an 
“S-curve”.  Figure  17  shows  two  such  calibration  curves,  one  for  a  point  source  and 
another  for  an  elliptical  Gaussian  beacon. 

Due  to  the  fact  that  the  curve  is  nonlinear,  the  sensor  is  typically  designed  to 
ensure  that  during  normal  operation  in  turbulence  most  of  the  local  tilt  on  a  sub- 
aperture  is  within  the  central,  linear  region  of  the  curve.  The  slope  of  the  curve  in 
this  region  determines  the  conversion  required  between  measured  and  applied  tilt. 
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In  this  research,  two  cases  are  compared.  The  initial  uses  the  traditional  method 
of  calibration,  with  a  diffraction  limited  spot  as  the  reference  source.  This  spot 
is  translated  across  the  detector  at  fine  increments,  and  the  centroid  calculation  is 
performed  at  each  location  to  determine  the  measured  tilt  and  produce  the  calibration 
curve.  The  gradient  of  the  linear  region  is  determined  by  using  MATLAB’s  polyf  it 
function  to  fit  a  first  order  polynomial  between  two  points  very  close  to  the  center  of 
the  curve.  The  second  case  performs  all  the  above  steps,  but  uses  the  more  realistic 
elongated  beacon  image  as  the  reference  source. 

The  relative  size  of  the  imaged  spot  and  the  pixels  that  measure  it  affect  the 
gradient  of  the  linear  region  of  the  calibration  curve.  It  is  therefore  during  sensor 
calibration  that  the  difference  between  calibration  with  a  traditional  reference  source 
and  an  elongated  source  will  be  highlighted,  and  provide  the  means  to  quantify  the 
effect  of  perspective  elongation  on  wavefront  measurements. 

3.3  Reconstruction 

To  reconstruct  the  aberrated  wavefront  that  is  represented  by  the  phase  screen 
described  in  Section  3.1.2,  the  local  tilt  across  each  subaperture  must  be  calculated 
and  applied  as  a  lateral  shift  of  the  Gaussian  function.  The  x  and  y  gradients  of  the 
high  resolution  phase  screen  is  determined  using  MATLAB’s  gradient  function.  The 
gradient  values  of  the  phase  screen  are  averaged  across  each  subaperture  to  determine 
the  average  wavefront  gradient  across  each  subaperture.  The  centroid  of  the  shifted 
images,  which  are  the  intensity-weighted  phase  gradients,  are  then  determined  and 
the  sensor  calibration  factor  is  applied. 

Once  slope  measurements  are  obtained,  the  original  wavefront  can  be  recon¬ 
structed  as  described  in  Section  2.1.3.  Typically,  calculating  the  reconstruction  matrix 
which  relates  the  measured  slopes  at  each  subaperture  to  the  measured  wavefront  at 
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32  subapertures  across  aperture 


Illuminated  apertures  are  active 

Figure  18.  WaveProp’s  subapdlg  toolbox,  used  to  generate  reconstruction  matrices  for 
different  system  geometries. 

each  of  the  DM  actuators  is  not  an  easy  endeavor.  Fortunately,  AOTools  contains  a 
tool  which  generates  the  geometry  matrix  and  various  reconstruction  matrices  for  a 
user-specified  AO  geometry.  This  tool  is  called  subapdlg,  and  it  allows  the  user  to 
input  the  number  of  actuators,  number  of  subapertures,  and  the  shape  of  the  pupil 
mask.  A  plot  of  the  geometry  used  in  this  research  is  shown  in  Fig.  18.  This  plot 
was  generated  by  subapdlg. 

Reconstruction  of  the  incident  wavefront  provides  enough  information  to  deter¬ 
mine  the  benefit  resulting  from  an  enhanced  calibration  method.  The  reconstructed 
wavefront  is  calculated  in  units  of  radians,  representing  the  phase  difference  between 
two  points  in  space.  The  wavefronts  of  the  original  and  enhanced  calibration  wave- 
fronts  can  be  compared  against  a  wavefront  that  is  reconstructed  from  the  original 
phase  screen  values  to  determine  the  wavefront  error  in  each  case.  The  root  mean 
square  error  of  the  phase  differences  at  each  subaperture  is  calculated  and  the  average 
over  all  active  subapertures  provides  a  metric  to  compare  the  methods. 
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IV.  Results  and  Analysis 


This  chapter  presents  the  results  that  were  obtained  using  the  methods  discussed 
in  Chapter  3.  The  impact  of  telescope  orientation,  turbulence  parameters,  and  de¬ 
tector  resolution  are  considered  to  characterize  the  effect  of  beacon  elongation  and 
rotation.  The  primary  metric  used  to  score  the  performance  of  calibration  methods 
was  Root  Mean  Square  (RMS)  phase  error  between  reconstructed  wavefronts. 

4.1  Rotation  Effects 

Equation  (18)  describes  how  to  determine  the  amount  of  rotation  caused  by  a 
Coude  path  for  a  given  trajectory.  Using  the  ISS  satellite  track  data  from  [18],  the 
Coude  rotation  angle  at  each  satellite  position  was  determined,  and  is  displayed  in 
Fig.  19.  For  this  particular  path,  the  telescope  aperture  rotates  at  21.4  deg/minute 
on  average,  and  rotates  less  at  lower  elevations  where  the  azimuthal  angle  to  the  ISS 
changes  much  less. 

4.2  Calibration  Effects 

As  described  in  Sec.  3.2.2,  the  wavefront  sensor  requires  calibration  to  determine 
the  detected  subaperture  tilt  that  is  caused  by  a  known  reference  tilt.  The  typical 
method  performs  the  calibration  using  a  point  source  as  the  reference  beacon,  which 
does  not  consider  the  effects  of  beacon  rotation  and  elongation.  Equation  (7)  describes 
that  the  resulting  PSF  using  square  lenslets  is  a  sine 2  function.  Figure  20  shows  the 
result  when  this  PSF  was  implemented  using  a  side  resolution  of  64  grid  points. 

The  width  of  the  PSF  is  an  important  factor  to  consider.  Since  the  literature  did 
not  describe  the  complete  specifications  of  the  sensor,  the  PSF  was  chosen  to  cover 
roughly  half  of  the  detector  area,  which  appears  to  match  empirical  data  provided 
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Figure  19.  Coude  rotation  caused  by  changes  in  target  azimuth  and  elevation. 
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Figure  20.  Diffraction  limited  point  spread  function.  This  axes  of  this  image  represent 
the  width  of  one  subaperture. 
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in  Ref.  [4],  The  validity  of  this  PSF  is  also  checked  against  metrics  used  by  tOSC 
to  design  Shack-Hartmann  sensors.  tOSC  define  a  parameter  ni0d,  or  ‘number  of 
spot  half- widths  per  subaperture’,  to  characterize  the  subaperture  field  of  view,  and 
defined  it  as 


d  d2 
riled  =  jx  =  jy 


(24) 


As  a  rule  of  thumb,  the  sensor  should  have  between  6  and  8  xni0d  across  each 
detector [7].  By  setting  the  null-to-null  width  of  the  diffraction  limited  PSF  in  the 
model  to  2  (in  arbitrary  units  due  to  the  undefined  lenslet  diameter  d  and  focal  length 
/),  the  diffraction  limited  spot  size  is  given  by [12] 


2/A 

d 

■ 

'  d 


2, 

1. 


(25) 


If  the  detector  width  is  set  to  3.5  (in  the  same  arbitrary  units),  n^,  is  also  equal 
to  3.5.  This  is  less  than  the  ideal  case  of  6-8  [7],  but  not  impractical.  Furthermore,  the 
number  of  pixels,  nptx  per  X/d  is  another  sensor  specification  relating  to  its  resolution. 
This  value  changes,  depending  on  the  detector  chosen,  but  typically  should  be  greater 
than  two  to  resolve  a  spot.  For  a  quad  cell  this  is  1.75,  and  for  a  4  x  4  detector  it  is 
3.5.  What  these  numbers  indicate  is  that  the  PSF  is  not  resolved  or  Nyquist  sampled 
in  a  quad  cell  detector,  but  is  for  the  4x4  scenario. 

Also,  a  PSF  at  this  width  is  generally  slightly  larger  than  the  geometric  beacon 
image  for  subapertures  near  the  laser  launch,  and  noticeably  smaller  than  the  geomet¬ 
ric  beacon  image  at  greater  elongations.  Figure  21  shows  the  minimum  and  maximum 
size  of  a  beacon  image  caused  by  perspective  elongation.  These  images  contain  the 
entire  FOV  of  the  subaperture,  as  is  the  case  of  Fig.  20,  for  easy  comparison.  The 
relative  sizes  of  the  beacon  image  and  PSF  determines  the  size  and  shape  of  the 
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Figure  21.  Beacon  images  at  two  subapertures.  The  left  image  shows  the  beacon  as 
viewed  from  the  subaperture  closest  to  the  laser,  and  the  right  image  shows  beacon 
from  the  furthest  subaperture.  These  beacon  images  are  taken  at  an  elevation  angle  of 
60  degrees,  and  would  be  slightly  different  at  other  angles,  as  discussed  in  Sec.  3.1.4 


diffraction  image  after  the  convolution  operation. 

4.2.1  Standard  Calibration. 

In  the  calibration  method  that  is  typically  performed,  the  detector  is  calibrated  by 
computing  the  centroid  of  the  PSF  as  it  is  translated  at  regular  intervals,  simulating 
the  effect  of  local  tilt.  The  calibration  curve  that  is  obtained  for  a  quad  cell  detector  is 
shown  in  Fig.  17.  The  y- axis  of  Fig.  17  is  measured  in  detector  width,  and  the  x-axis 
in  radians  of  wavefront  slope,  which  can  be  easily  converted  to  optical  path  difference 
measured  in  wavelengths.  Note  that  the  maximum  measurable  displacement  for  a 
quad  cell  detector  is  ±0.25  x  d,  where  d  is  the  detector  width.  This  is  the  center  of 
the  furthest  pixel,  corresponding  to  the  case  where  the  entire  PSF  is  contained  within 
this  pixel.  At  this  tilt,  the  sensor  is  beyond  its  useful  linear  region,  and  the  sensor 
is  designed  such  that  the  expected  magnitude  of  wavefront  slopes  do  not  cause  this 
amount  of  tilt. 

This  calibration  curve  also  determines  the  extent  of  the  linear  region  of  the  detec- 


39 


tor.  In  Fig.  17,  this  corresponds  to  roughly  30%  of  the  detector  width.  The  sensor 
is  designed  such  that  the  maximum  expected  local  tilt  rarely  deflects  the  centroid  of 
the  PSF  beyond  this  point.  This  ensures  that  a  linear  calibration  factor  can  be  used 
to  represent  the  curve  without  too  much  error. 

To  determine  the  calibration  factor,  MATLAB@’s  polyfit  command  was  used 
to  fit  a  line  to  the  center  region  of  the  curve  as  shown  in  Fig.  17.  Due  to  symmetry 
of  the  point  source’s  image,  this  curve  is  identical  for  both  the  x  and  y  directions. 

4.2.2  Enhanced  Calibration. 

When  considering  the  effects  of  elongation  and  rotation  on  the  beacon  image 
shape,  a  new  reference  source  is  required  for  each  subaperture.  This  alters  the  cali¬ 
bration  curve  at  each  subaperture,  since  the  beacon  image  varies  depending  on  the 
relative  distance  and  angle  between  the  subaperture  and  the  laser.  Creating  a  refer¬ 
ence  source  that  emulates  beacon  elongation  and  rotation  for  laboratory  use  is  difficult 
and  has  not  currently  been  achieved.  This  research  explores  a  low  risk  solution  to 
providing  these  calibration  values  without  producing  a  new  reference  source. 

The  slope  of  the  line  fitted  to  the  linear  region  of  the  calibration  curves  (the  cal¬ 
ibration  factor)  for  every  subaperture  is  shown  in  Fig.  22.  The  lower  calibration 
factors  represent  subapertures  where  the  beacon  has  greater  elongation,  which  gen¬ 
erally  increases  the  linear  region  of  the  detector.  Additionally,  the  calibration  factors 
for  the  x  and  y  directions  differ  due  to  the  anisotropic  beacon  image. 

These  new  calibration  factor  predictions  are  the  key  to  reducing  the  phase  error 
of  the  reconstructed  wavefront.  If  the  gradients  from  a  standard  calibration  (which 
does  not  consider  beacon  elongation)  were  displayed  in  a  similar  manner,  it  would 
remain  a  constant  value  across  the  detector  due  to  symmetry  of  the  PSF.  Furthermore, 
the  Coude  rotation  at  the  sensor  causes  these  calibration  gradients  to  change  with 
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Figure  22.  Calibration  factors  over  the  entire  sensor.  Note  the  relative  direction  of 
the  laser  and  its  effect  on  the  calibration  factor  across  the  aperture.  Subapertures 
furthest  from  the  laser  are  affected  by  the  most  elongation,  lowering  the  gradient  of 
the  calibration  curve. 


telescope  orientation,  preventing  the  use  of  a  simple  calibration  method  using  an 
elongated  source  in  one  orientation. 

4.3  Sensing  Turbulence 

To  model  the  effect  of  turbulence,  randomly  drawn  phase  screens  are  generated 
using  the  AOTools  function  kolmogphzscreen,  such  as  the  phase  screen  shown  in 
Fig.  23.  This  phase  screen  is  continually  used  throughout  this  chapter  to  illustrate 
the  benefit  of  the  enhanced  calibration  method. 

To  generate  this  phase  screen,  1024  grid  points  were  used  across  an  aperture  size 
of  3.5m.  The  coherence  diameter  is  set  to  6.1cm  for  the  standard  case,  and  50%  and 
200%  of  this  value  for  comparison,  based  upon  r0  calculations  for  the  SOR  site,  as 
described  in  Sec.  3.1.2.  The  phase  screen  is  generated  without  adding  any  global  tilt, 
similar  to  a  wavefront  after  the  majority  of  tilt  has  been  removed  by  a  fast  steering 
mirror. 

The  local  x  and  y  wavefront  gradients  are  then  computed  using  MATLABlgj’s 
gradient  function.  These  gradients  are  subaperture-averaged  to  align  to  WFS  mea¬ 
surements,  as  shown  in  Fig.  24.  These  local  gradients  are  converted  to  a  translation 
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Figure  23.  Phase  screen  simulating  the  effects  of  turbulence  on  a  planar  wavefront. 


of  the  PSF  at  each  subaperture,  modeling  the  effect  of  propagation  through  a  Shack- 
Hartmann  lenslet  to  the  detector.  The  shift  is  induced  as  the  x0ffset  and  y  off  set 
variables  in  Eq.  (22).  The  displacement  of  the  turbulent  Shack-Hartmann  spots  are 
then  measured  using  the  centroid  algorithm  over  the  desired  number  of  pixels. 

The  wavefront  slope  at  each  subaperture  is  then  determined  by  applying  the 
centroid-dependent  calibration  factor  to  each  centroid  measurement.  As  expected, 
the  measured  slopes  obtained  from  the  standard  calibration  varied  from  the  enhanced 
calibration. 

4.4  Reconstruction 

Finally,  the  measured  wavefront  slopes  are  reconstructed  into  phase  values  through 
Eq.  (17).  The  reconstructor  matrix,  H,  was  obtained  through  AOTools  subapdlg 
tool.  Several  stages  of  the  subaperture  averaging  and  reconstruction  process  are  shown 
in  Fig.  25.  Although  both  results  appear  similar  to  the  original  phase  measurements, 
taking  the  difference  between  them  shows  more  clearly  how  much  more  accurate 
the  enhanced  calibration  method  is.  Figure  26  shows  the  difference  between  the 
reconstructed  and  original  phase.  There  is  a  fixed  color  scale  for  both  images  in 
Fig.  [26],  highlighting  the  much  lower  phase  error  produced  by  using  the  enhanced 
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Quiver  plot  of  wavefront  gradients 


Figure  24.  Quiver  plot  of  randomly  generated  wavefront  slopes  taken  from  the  phase 
screen  in  Fig.  23  before  any  calibration  factor  has  been  applied.  The  magnitude  of 
the  largest  slope  is  8.32  radians  per  subaperture.  These  arrows  also  correspond  to  the 
relative  displacement  of  Shack-Hartmann  spots. 


calibration  method. 

4.5  Performance 

The  metric  used  to  determine  the  performance  of  the  model  is  the  root  mean 
square  error  (RMSE)  between  the  incident  and  reconstructed  wavefront  phase.  To 
eliminate  the  effect  of  reconstruction  error,  the  reconstructed  phases  are  not  com¬ 
pared  directly  to  the  original  phase  screen,  but  instead  to  the  reconstruction  of  the 
downsampled  phase  screen  gradients.  50  iterations  using  different  atmospheric  phase 
screens  were  used  at  each  telescope  orientation  corresponding  to  a  satellite  position, 
and  the  RMS  phase  error  at  each  subaperture  was  averaged  over  each  iteration  to 
determine  the  final  RMSE. 

The  number  of  iterations  ran  for  each  orientation  of  the  telescope  was  verified 
by  calculating  the  standard  deviation,  (Jrmse,  and  mean,  Urmse  RMS  phase  error 
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Figure  25.  Example  wavefront  at  various  stages  of  reconstruction.  Initially,  the  wave- 
front  is  downsampled  from  a  high  resolution  to  be  fitted  to  the  DM.  The  difference 
between  reconstructed  phase  from  the  different  calibration  techniques  is  only  slightly 
evident  here. 
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Figure  26.  Difference  between  reconstructed  and  original  phase. 
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Figure  27.  Checking  that  enough  iterations  of  turbulence  are  used  in  this  model. 


Table  1.  Average  Root  Mean  Square  Error  of  reconstructed  phase  over  ISS  path. 


d/r0: 

3.59 

1.79 

0.90 

detector  size: 

2x2 

4x4 

2x2 

4x4 

2x2 

4x4 

standard  cal 

9.59 

2.94 

3.89 

0.63 

1.63 

0.21 

enhanced  cal 

7.17 

2.75 

1.98 

0.50 

0.29 

0.13 

error  reduction 

25.2% 

6.5% 

49.1% 

20.6% 

82.2% 

38.1% 

after  each  iteration.  The  results  of  this  analysis  with  a  quad  cell  detector  with  using 
standard  calibration  and  with  r o  =  3.05cm  are  shown  in  Fig.  27.  As  a  general  rule, 
the  number  of  iterations  was  deemed  sufficient  when  the  standard  deviation  of  the 
data  set  was  less  than  10%  of  the  mean.  In  this  circumstance,  50  iterations  is  more 
than  sufficient.  Other  configurations  were  tested  with  similar  results. 

The  average  phase  RMSE  across  all  detectors  over  50  iterations  of  the  ISS  path  is 
shown  in  Table  1,  as  a  function  of  r o  and  detector  resolution.  As  expected,  the  results 
show  a  decrease  in  RMS  phase  error  as  the  d/r0  value  decreases,  since  the  system  is 
affected  by  weaker  turbulence.  Similarly,  the  transition  from  a  quad  cell  detector  to  a 
4x4  pixel  detector  reduces  the  RMS  phase  error  in  each  case.  The  RMSE  is  reduced 
by  49.1%  through  using  the  enhanced  calibration  method,  compared  to  20.6%  if  the 
detector  was  upgraded  from  a  quad  cell  to  a  4  x  4  pixel  detector.  Detectors  with  even 
greater  resolution  were  investigated,  but  the  reduction  in  error  was  almost  negligible. 
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Figure  28.  Effect  of  turbulence  on  sensor  measurements. 


4.6  Analysis 

The  quad  cell  detector  with  r0  =  6.1cm  was  considered  the  baseline  configuration 
against  which  other  configurations  were  compared.  The  value  of  r0  was  first  increased 
then  later  decreased  by  a  factor  of  2  to  determine  the  effect  of  atmospheric  turbulence 
strength  on  the  results.  Secondly,  the  detector  resolution  was  improved,  and  finally, 
a  second  satellite  path  was  used  to  validate  results.  Overall,  the  results  show  that 
considering  elongation  effects  during  calibration  can  reduce  RMS  wavefront  error 
by  approximately  50%  for  the  standard  operating  configuration,  with  typically  less 
noticeable  benefits  outside  of  this  configuration.  Images  of  the  complete  results  are 
displayed  in  Sec.  4.7.3. 

4.7  Atmospheric  Effects 

The  atmospheric  coherence  diameter  was  varied  to  represent  varying  turbulence 
strengths  to  determine  the  effect  on  wavefront  error.  It  was  expected  that  the  model 
would  show  improved  performance  with  larger  coherence  diameters.  The  mean  RMSE 
of  the  reconstructed  wavefront  is  shown  in  Fig.  28. 

As  expected,  increasing  the  atmospheric  coherence  diameter  decreased  RMS  phase 
error.  In  the  quad  cell  configuration,  RMSE  was  reduced  from  from  7.17  rad  at  d/r^  = 


46 


2x2  cell,  standard  cal 


10 

20 

30 


10  20  30 

RMSE  =  3.8926 


2x2  cell,  enhanced  cal 


10  20  30 

RMSE  =  1.9847 


4x4  cell,  standard  cal  4x4  cell,  enhanced  cal 


10  20  30  10  20  30 

RMSE  =  0.63376  RMSE  =  0.50091 


Figure  29.  Effect  of  detector  resolution  on  sensor  measurements.  The  left  column  show 
the  effect  of  improved  resolution  without  using  the  enhanced  calibration  method,  while 
the  right  column  shows  the  improvement  after  applying  the  appropriate  calibration 
factors. 


3.59  to  0.29  rad  at  d/ro  =  0.90.  This  is  a  positive  result,  confirming  expectations. 

4.7.1  Detector  Resolution. 

The  detector  was  increased  in  resolution  from  a  quad-cell  to  a  4  x  4  cell  detector 
to  characterize  the  effects  of  detector  resolution  on  measurements.  The  improvement 
gained  by  using  the  enhanced  calibration  technique  was  noticeably  less  in  this  case 
compared  to  the  quad  cell  detector,  due  to  the  extended  linear  region  for  both  elon¬ 
gated  and  diffraction  limited  spots  over  4  pixels.  Detectors  with  resolutions  greater 
than  4x4  pixels  were  briefly  investigated,  but  showed  almost  negligible  improvement 
due  to  a  much  greater  linear  region  and  much  more  precise  centroid  calculations. 
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Table  2.  Average  Root  Mean  Square  Error  of  reconstructed  phase  over  GPS  path. 


d/ro. 

3.59 

1.79 

0.90 

detector  size: 

2x2 

4x4 

2x2 

4x4 

2x2 

4x4 

standard  cal 

9.52 

3.05 

3.91 

0.66 

1.63 

0.23 

enhanced  cal 

7.01 

2.83 

1.94 

0.51 

0.27 

0.14 

error  reduction 

26.4% 

7.2% 

50.4% 

22.7% 

83.4% 

39.1% 

4.7.2  Trajectory  Differences. 

The  GPS  satellite  orbit  path  showed  little  difference  in  the  effect  on  performance 
of  the  sensor  when  compared  to  that  of  the  ISS.  This  was  expected,  as  the  overall 
sensor  performance  should  not  be  affected  by  changes  in  the  Coude  rotation  of  the 
beacon.  The  full  results  of  the  GPS  trajectory  are  shown  in  Table  2. 

4.7.3  Complete  Results. 

The  complete  results  for  each  turbulence  strength,  detector  resolution  and  path 
type  combination  are  shown  in  Figs.  30-33.  These  figures  show  the  difference  in 
RMSE  between  the  typical  calibration  method  (top)  and  the  enhanced  calibration 
method  (bottom).  The  d/ro  value  decreases  from  left  to  right.  For  a  3.5m  telescope 
with  32  subapertures  across,  d  =  10.94cm.  The  r0  value  for  the  center  plots  is  6.1cm, 
but  is  50%  smaller  for  the  left  plots  and  100%  greater  for  the  right.  These  figures 
show  that  the  enhanced  calibration  method  produces  less  error  than  the  standard 
calibration  in  all  scenarios,  but  provides  the  greatest  benefit  in  the  quad  cell  scenario. 
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Figure  30.  Full  results,  ISS  path,  quad  cell. 
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Figure  31.  Full  results,  ISS  path,  4x4  cell. 
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Figure  32.  Full  results,  GPS  path,  quad  cell. 
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Figure  33.  Full  results,  GPS  path,  4x4  cell. 
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V.  Conclusions  and  Future  Work 


This  research  has  developed  a  model  to  represent  the  current  AO  system  at  the 
Starfire  Optical  Range.  Although  this  model  did  not  consider  every  aspect  that 
affects  the  AO  system,  including  all  hardware  specifications  and  noise  sources,  it  still 
demonstrated  the  improvement  in  measurement  capability  that  can  be  achieved  by 
considering  Coude  rotation  and  perspective  elongation. 

This  chapter  outlines  some  of  the  areas  that  would  require  improvement  before  this 
model  could  be  implemented  into  a  current  system,  and  also  some  of  the  interesting 
features  that  were  discovered  through  the  development  of  this  model. 

5.1  Desired  Improvements 

The  model  that  was  developed  was  kept  flexible  in  its  design  to  easily  compare 
different  scenarios.  Some  areas  of  its  design  were  kept  generic  for  simplicity,  and 
require  alteration  before  applying  the  model  to  a  specific  system.  For  example,  the 
PSF  of  the  real  optical  equipment  should  be  used  to  ensure  that  it  correctly  matches 
the  system.  Similarly,  the  reconstruction  matrix  used  to  convert  the  wavefront  slope 
measurements  into  DM  commands  was  developed  using  AOTools  and  fits  the  data 
onto  a  33  x  33  actuator  DM.  The  actual  reconstruction  matrix  for  the  real  DM  should 
be  used  in  its  place.  Also,  the  AO  system  components  were  modeled  as  being  ideal, 
and  there  was  no  accounting  for  the  spatial  influence  function  clue  to  neighboring 
actuators  on  the  DM. 

The  resolution  of  the  beacon  images  at  each  subaperture  was  chosen  to  be  64  x 
64  grid  points,  as  a  compromise  between  the  time  it  took  to  run  each  simulation 
and  the  quality  of  the  results.  This  was  chosen  due  to  the  large  number  of  test 
cases  that  were  analyzed,  but  this  resolution  should  be  improved  to  achieve  more 
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accurate  results  once  the  system  has  been  modeled  correctly.  Additionally,  the  beacon 
images  themselves  are  merely  approximations,  and  further  analysis  into  the  effects  of 
perspective  elongation,  such  as  experimental  data,  is  essential  to  provide  a  closer  fit 
to  real  effects.  Other  techniques  such  as  image  correlation  may  be  a  suitable  option 
to  account  for  the  elongation. 

Finally,  implementing  this  model  in  a  laboratory  would  be  of  tremendous  benefit 
in  validating  the  model.  Unfortunately,  realistic  laboratory  sized  extended  beacons 
are  not  currently  available. 

5.2  Implementation 

If  this  model  was  to  be  implemented,  it  is  unlikely  that  improved  calibration 
factors  could  be  determined  in  real-time  using  current  computer  processing  capabili¬ 
ties.  To  overcome  this,  two  methods  could  be  implemented.  The  first  would  require 
pre-calculating  calibration  factors  using  the  predicted  trajectory  of  the  object  to  be 
viewed.  This  would  be  relatively  simple  to  implement,  but  would  require  a  new  cali¬ 
bration  model  for  each  trajectory.  This  could  be  done  either  between  satellite  passes 
if  there  is  sufficient  time  or  earlier  in  the  day  once  the  nightly  target  list  has  been 
determined. 

An  alternate  solution  could  be  to  run  the  model  for  fixed  positions  in  azimuth 
and  elevation  to  cover  one  quarter  of  the  visible  sky,  as  shown  in  Fig.  34.  The 
remaining  sky  would  be  determined  by  rotating  the  x  and  y  axes  of  the  original 
measurement.  Provided  a  sufficient  number  of  points  are  analyzed  in  the  space,  the 
remaining  calibration  factors  could  be  interpolated  from  nearby  measurements.  It  is 
expected  that  these  calibration  factors  could  be  accessed  in  real-time  scenarios,  and 
would  only  require  measurements  during  the  initial  implementation  of  this  technique. 
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Figure  34.  Quadrant  of  the  sky  showing  an  example  of  how  calibration  could  be  deter¬ 
mined  at  specific  azimuth-elevation  angles  and  interpolated  for  angles  between  these 
points. 

5.3  Key  Results 

This  research  effectively  showed  that  using  standard  calibration  techniques  in  the 
presence  of  elongated  references  sources  introduces  inaccuracies  in  wavefront  mea¬ 
surement.  The  model  developed  through  this  research  predicts  that  errors  in  the 
reconstructed  phase  can  be  reduced  by  up  to  50%  by  correctly  predicting  the  shape 
of  the  measured  Shack-Hartmann  spot. 

Factors  such  as  d/r0  were  characterized,  revealing  that  smaller  values  of  d/r0 
produced  less  errors  in  the  reconstructed  phase,  which  is  expected.  This  was  a  useful 
technique  to  verify  the  model.  Additionally,  the  effect  of  detector  resolution  was 
investigated,  and  again  the  results  matched  expectations.  This  research  determined 
that  improving  the  detector  from  a  quad  cell  to  a  4  x  4  pixel  detector  improved  results 
more  than  by  simply  accounting  for  the  source  elongation  during  sensor  calibration. 
It  is  expected  that  this  greater  improvement  will  come  with  a  much  greater  financial 
cost,  in  addition  to  the  complexity  associated  with  its  implementation. 

One  surprising  result  is  the  location  of  phase  error.  Despite  the  gradient  mea¬ 
surements  visibly  being  most  inaccurate  at  subapertures  affected  by  the  greatest  per- 
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spective  elongation,  this  did  not  relate  to  the  greatest  phase  errors  being  present  at 
these  locations  after  phase  reconstruction.  Analysis  of  the  location  of  the  phase  error 
relevant  to  distance  from  the  laser  indicated  that  there  was  little  to  no  correlation 
between  phase  measurement  error  and  distance  to  laser. 

One  possible  explanation  for  this  is  that  the  results  were  obtained  using  a  tilt- 
removed  reconstructor,  whereas  there  may  still  be  some  residual  tilt  present  in  the 
wavefront,  even  after  being  corrected  by  a  FSM.  This  removal  of  tilt  may  mask  the 
location  of  errors. 

5.4  Future  Work 

Further  work  on  this  model  is  highly  encouraged,  and  recommendations  on  how 
to  improve  the  model  are  provided  in  this  section.  These  include  incorporating  the 
enhanced  calibration  factors  into  the  reconstruction  matrix,  replacing  the  beacon  im¬ 
age  with  a  more  realistic  image,  and  adjusting  atmospheric  parameters  with  telescope 
orientation. 

5.4.1  Dynamic  Reconstructor. 

In  the  current  implementation  of  the  model,  measured  wavefront  slopes  are  mul¬ 
tiplied  by  their  respective  calibration  factors  and  then  reconstructed  using  Eq.  (17). 
A  slight  variation  of  this  method  would  be  to  post-multiply  the  reconstructor  matrix 
H  with  the  matrix  form  of  the  calibration  factors,  C,  such  that 

0  =  HCs.  (26) 

Despite  achieving  the  same  results,  performing  this  operation  would  allow  analysis  of 
the  effect  of  elongation  and  pupil  rotation  on  the  new  reconstruction  matrix  HC,  and 
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could  possibly  be  a  precursor  to  developing  a  predictive  and  dynamic  reconstruction 
matrix. 

5.4.2  Improved  Beacon  Image. 

Another  improvement  that  could  be  readily  implemented  involves  replacing  the 
first-order  approximation  of  the  beacon  image  described  by  Eq.  (22)  with  a  more  real¬ 
istic  image.  The  beacon  image  could  be  created  by  considering  the  three-dimensional 
nature  of  the  sodium  layer,  including  the  varying  diamter  of  the  laser  beam  as  it 
propagates  and  the  Gaussian  distribution  of  sodium  atoms  within  the  layer.  The 
image  could  be  modeled  using  multiple  two-dimensional  slices  through  the  sodium 
layer,  and  the  resulting  beacon  propagated  to  the  telescope  using  wave  optics  soft¬ 
ware.  Furthermore,  the  calibration  method  using  the  Gaussian  approximation  could 
remain  unchanged  to  characterize  the  sensitivity  of  the  results  in  this  research  to 
changes  in  the  shape  of  the  beacon. 

5.4.3  Varying  Atmospheric  Parameters. 

In  this  model,  the  same  C\  profile  and  therefore  ro  value  was  used  to  model 
atmospheric  turbulence,  regardless  of  the  elevation  angle  of  the  telescope.  This  was 
done  to  characterize  sensor  performance  at  varying  elevation  angles  due  to  beacon 
elongation  and  rotation  separately  from  the  effects  of  varying  atmospheric  parameters. 
A  more  realistic  model  would  adjust  the  turbulence  parameters  with  the  elevation 
angle  of  the  telescope. 

5.5  Impact  of  Work 

The  cumulative  effect  of  pupil  rotation  and  perspective  elongation  on  wavefront 
sensor  measurements  is  a  unique  issue  affecting  the  Starfire  Optical  Range.  Despite 
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currently  performing  well,  the  system  could  still  be  improved  by  simply  implementing 
the  methods  outlined  in  this  thesis.  Although  once  more  complex  factors  are  consid¬ 
ered,  the  actual  reduction  in  RMSE  of  the  measured  phase  may  not  be  as  great  as 
determined  by  this  model  («50%),  a  definite  improvement  is  still  achievable. 
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Appendix  A.  MATLAB  Code 


Listing  A.l.  performance. m 


1  function  [EnhCal  StdCal]  =  performance (nloops,  rOfac,  PIXres,  az,  el) 

2  %PERFORMANCE  determines  the  measurement  error  for  standard  and 

3  %enhanced  calibration  methods  with  LGS  elongation. 

4  %  Designed  to  run  several  test  cases  with  various  rO  values  and 

5  %detector  resoltions.  Produces  EnhCal  and  StdCal,  the  error  using 

6  %each  random  phase  screen  compared  against  an  ideal  system. 

7  %NOTE :  Results  are  straight  differences  between  ideal  reconstructed 

8  %phase  and  the  two  calibration  methods.  Further  analysis  such  as 

9  %taking  the  root  mean  square,  and/or  averaging  over  the  active 

10  %subapertures  (877  for  32x32)  must  still  be  performed! 

11  %NOTE :  Don't  forget  to  generate  H,  the  reconstruction  matrix  from  the 

12  %subapdlg  toolbox!! 

13  % 


14  %  @author: 

FLTLT  Russell  McGuigan 

is  %  (Mate: 

21Feb2011 

16  %  @ inputs: 

17  %  nloops 

number  of  random  phase  screen  draws 

is  %  rOfac: 

factor  to  alter  rO  (for  comparisons) 

19  %  PIXres 

detector  resolution  (e.g.  2  for  2x2  quad  cell) 

20  %  a  z  /  e  1 : 

azimuth/elevation  angles  [rad] 

21  %  Mutputs 

22  %  EnhCal 

error  of  enhanced  calibration  method 

23  %  StdCal 

error  of  standards  calibration  method 

2 4  ooooooooooooooooooooooooooooooooooooooooo 

25 

2o  oooooooooooooooooooooo 

27  %  Initial  Setup 
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Zo  oooooooooooooooooooooo 


29 

Diameter= 

3.5; 

%Telescope  diameter,  [m] 

30 

Nsub= 

32; 

%Number  of  subapertures  across 

telescope 

31 

L_sep= 

i; 

%radial  separation  of  laser  from  edge  of  telescope, [m] 

32 

theta_0= 

0; 

%initial  offset  of  Coude  path 

[deg] 

33 

34 

r0=0 . 061*r0fac; 

%Fried's  parameter 

35 

N=2"10; 

%Resolution  of  phase  screen 

36 

PSFres=64 ; 

%resolution  of  PSF  and  Ellipse 

image 

37 


38  myPath  =  myTrack (az, el) ;  %generate  track  data 


39  myLayer  =  myLayers (myPath) ; %determine  distance  to  sodium  layer 


40  %boundaries  for  given  track  data 


4i  myTelescope  =  myTel (  Nsub, Diameter , L_sep , [],  PIXres,  theta_0); 


42  %generates  subaperture  grid 


43  myLaser=myC_rotate (myTelescope, myPath) ; 


44  %Determines  Coude  rotation  for  a  given 


45  %path  and  telescope 


46  myTurb  =  myAtmosphere (Nsub, Diameter , rO , N) ; 


47  %Generates  random  phase  screen 


48  PSFwidth=l;  %Determines  width  of  sinc~2  function 


49  [PSF  eta  nu] =myPSFgen (myTelescope, PSFres, PSFwidth) ; 


so  %Generates  PSF  for  a  given  system 


5i  minor_ax=0 . 1 ;  %Width  of  ellipse  minor  axis 


52  res=8;  %factor  to  upsample  calibration  data 


53  %to  achieve  smooth  curves 


54  shiftfac=0 . 09*PSFres/32;  %Factor  to  translate  images  due  to 


55  %local  tilt.  Similar  effect  to 


56  %changing  the  speed  of  a  lens 
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oo  oooooooooooooooooooooooooooooo 


59  %%  Standard  Calibration 
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60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
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82 

83 

84 

85 

86 

87 

88 

89 
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2'9'2'2-2'9-9'2'9'2'9'2'9'2'9-9'9-9'2-9'2'9'2'9-2'2-9'9-9'2' 

oooooooooooooooooooooooooooooo 

calpos= [— 15  :  3 : 15 ] ;  %radian  range  to  calibrate  sensor  over 

myLayer2=myLayer ;  %reproduce  layer  values,  used  to 

myLayer2 .Mean=0;  %produce  a  'null'  ellipse  for  standard 

%calibration 

for  xcal=l : length (calpos ) ;  %determines  centroid  of  each  input 
ycal=0;  %translates  in  x  direction  only 

ellipse=myEllipseGen (my Laser , 1 ,  my Tele scope . x ( 1 , 1 ) , . . . 
myTelescope . y ( 1 , 1 ) , eta,  nu,  myLayer2, 0 . 005, 1, . . . 
calpos (xcal) , ycal, shiftfac)  ; 

%generate  ellipse  for  one  subaperture. 

%based  on  telescope/target  parameters 
convolved_ellipse=conv2 (PSF,  ellipse)  ; 

%convolves  ellipse  and  PSF 

cut of f=round ( 0 . 25*length ( convolved_ellipse ) ) ; 

%determines  excess  padding  to  remove 
conv_ellipse4=convolved_ellipse (cutoff : 3*cutoff,  . . . 

cutof f : 3*cutof f ) ;  %remove  padding  from  convolution 

integrated_convolved=myInt (conv_ellipse4 , myTelescope . res) ; 
%integrate  images  over  detector  pixels 
% (effectively  determines  photocounts) 

sum_tot=sum ( sum ( integrated.convolved) ) ;  %total  photocounts 
nocahx  (xcal)  =sum  (sum  (integrated_convolved.  *  .  .  . 

myTelescope  .  x_loc/sum_tot ) ) /myTelescope . Diameter; 
%determines  centroid  of  image 

end 

vec_x_2a=interp (nocahx, res) ; %interpolates  measurements  for  a 
%smooth  cal  curve 

vec_x_2=vec_x_2 a ( 1 : end— res+1 ) ; %removes  final  values  for  same  axis 
%plotting 

steps2=length ( vec_x_2 ) ;  %determines  length  of  cal  curve 

x_ax_2=linspace (calpos ( 1 ), calpos (end) , steps2 );  %used  to  aligns  cal 
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93 

94 
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%curve  to  input  value 

mid=round (length (vec_x_2 ) /2 ); %determines  centre  cal  point 
rng=round ( 0 . 07 ^length ( vec_x_2 ) ) ; %determines  range  over  which  to 
%fit  linear  region  (0.07=  +and— 7%) 

test_x=polyf it (x_ax_2 (mid— rng : mid+rng) ,  vec_x_2 (mid— rng : mid+rng) , 1) ; 
%fit  linear  polynomian  (line)  to  curve 

nocal_x2=test_x (1) ;  %takes  gradient  of  line  as  cal  factor 

nocal_y2=nocal_x2 ;  %x  and  y  are  symmetric 

2'9'2'2-2'9-9'2'9'2'9'2'9'2'9-9'9'9'2'9'2'9'2'9-2'2-9'9-9'2'9'2'9'2'2-2'9-2'9'9'2'9'2'9-2'9-2'2-9'9-9'2'9'2'9' 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%%  Enhanced  Calibration 

2'2-2'9'2'9-9'2'9'2'2-9'2-2'2'9'9'9'2'9'2'2-2'2'2'9'9'9-9'2'9'2'2-9'9'2'9'9'9-9'2'9'2'2'2'2'2-9'9'9-9'2'9'2'2- 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%Performs  same  calibration  as  standard  cal,  but  generates  a 
%unique  beacon  images  at  each  subaperture  to  model  actual 
%measurements 

for  x_pos  =  l:Nsub  %Perform  for  all  subapertures 

for  yrpos  =  l:Nsub  % (within  mask) 

if  (myTelescope .mask(x_pos,y_pos)>0) 
for  xcal=l : length (calpos) ; 
ycal=0 ; 

ellipse=myEllipseGen2 (myLaser, . . . 

1 ,  myTelescope .x(x_pos,y_pos) ,  .  .  . 
myTelescope . y (x_pos , y_pos ) , eta,  nu, . . . 
myLayer, minor_ax, 1, calpos (xcal)  ,  ycal,  .  .  . 
shiftfac) ; %generate  unique  ellipse 
convolved_ellipse=conv2 (PSF, ellipse) ; 
cutof f= (0 . 25*length (convolved.ellipse) ) ; 
conv_ellipse4=convolved_ellipse (ceil (cutoff) : . . . 
floor  (3*cutof f ) , ceil (cutoff)  : . .  . 
floor  (3*cutof f ) ) ; 

integrated_convolved=myInt ( conv_ellipse4 , . . , 
myTelescope . res ) ; 

sum_tot=sum ( sum ( integrated_convolved) ) ; 
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mycal.x  (xcal)  =sum  (sum  ( integratecLconvolved .  *  .  .  . 
myTelescope  .  x_loc/ sum_tot )  )  .  .  . 

/myTelescope . Diameter; 

end 

for  ycal=l : length (calpos )  %x  and  y  cals  will  differ 
xcal=0 ; 

ellipse=myEllipseGen2 (myLaser, . . . 

1  ,  myTelescope  .x(x_pos,y_pos)  ,  .  .  . 
myTelescope .y (x_pos, y_pos) , . . . 
eta,  nu,  myLayer, minor_ax,  1,  . . . 
xcal, calpos (ycal) , shiftfac) ; 
convolved_ellipse=conv2 (PSF, ellipse) ; 
cutoffs (0 . 25*length ( convolved_ellipse ) )  ; 
conv_ellipse4=convolved_ellipse  (ceil  (cutoff)  :  .  .  . 
floor (3*cutoff ) , ceil (cutoff) : . . . 
floor (3*cutoff ) ) ; 

integrated_convolved=myInt  ( conv_ellipse4 ,  .  .  . 
myTelescope . res ) ; 

sum_tot=sum  ( sum  ( integrated_convolved)  )  ; 
mycal.y  (ycal)  =sum  (sum  ( integrated.convolved .  *  .  .  . 
myTelescope  .  y.loc/ sum.tot )  )  .  .  . 

/myTelescope . Diameter; 

end 

vec_x_2a=interp (mycal_x,  res)  ; 

vec_x_2=vec_x_2a  ( 1 :  end— res+1 )  ; 

vec_y_2a=interp (mycal_y,  res)  ; 

vec_y_2=vec_y_2a  (1 :  end— res  +  1)  ; 

steps2=length ( vec_x_2  )  ; 

x_ax_2=linspace (calpos (1) , calpos (end) , steps2) ; 

mid=round ( length ( vec_x_2  )  / 2)  ; 

rng=round ( 0 . 1* length ( vec_x_2 ) )  ; 
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t est _x=polyf it (x_ax_2 (mid— rng :mid+rng) ,  . . . 

vec_x_2  (mid— rng : mid+rng) , 1) ; 
t est _y=polyf it (x_ax_2 (mid— rng :mid+rng) , . . . 

vec_y_2 (mid— rng : mid+rng) , 1) ; 
cal_x  (x_pos,y_pos)=test_x  (1)  ;%store  cal  factor 
cal_y  (x_pos  ,  y_pos  )  =test_y  (1)  ; 

else 

cal_x  (x_pos  ,  y_pos  )  =0; 
cal_y  (x_pos  ,  y_pos  )  =0; 

end 

end 

end 

2'9'9'9-2'9-9'2'9'2'9'2'9'2'9-9'9-9'2'9'2'9'2'9-2'9-9'9-9'2-9'2'9'2'9-2'9-9'9-9'2' 

ooooooooooooooooooooooooooooooooooooooooo 

%%  Apply  a  random  phase  screen  to  score  performance 

2'9'2'9-2'9'9'2'9'2'9'2'9'2'9-9'9-9'2'9'2'9'2'9-2'9-9'9-9'2-9'2'9'2'9-2'9-9'9-9'2' 

ooooooooooooooooooooooooooooooooooooooooo 

load  H  %load  reconstructor  matrix  from  subapdlg  toolbox 
%must  have  pre  saved  and  contain  ActiveSubs  and 
%Active  acts.  Run  toolbox  and  use  command: 

% 1  save  H  H  HOtlt  ActiveActs  ActiveSubs’ 

H=HOtlt;%Use  tilt  removed  reconstructor 

EnhCal=zeros ( [nloops, Nsub+1 , Nsub+1 ] ) ; %pre— allocate  result  arrays 
StdCal=zeros ( [nloops, Nsub+1 , Nsub+1 ] ) ; 
for  idxprime=l : nloops  %perform  for  each  iteration 
phz=myTurb .phase;  %extract  phase  screen  data 
phz_ds=interp2 (phz ,  1 inspace ( 1 ,  N,  Nsub+1 )  ,  .  .  . 
transpose ( 1 inspace ( 1 , N, Nsub+1 ) ) ) ; 

%downsample  phase  screen  for  display,  if  required 
steps=round ( 1 in space ( 1 , length (phz ) , Nsub+1 ) ) ; 

screen_grad_x=myTurb . screen_grad_x;  %extract  gradient  values 
screen_grad_y=myTurb  .  screen_grad_y; 

mask=logical (ActiveActs ); %converts  matrix  to  logicals 
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S=  [  screen_grad_y  (Active Subs)  ;  scr een.gr acLx  (Active Subs)  ]  ; 
%creates  ideal  slope  matris 

P=zeros ( size (Act iveActs )) ;  %initialise  reconstructed  phase 
P(mask)=H*S;  %Rreconstruct s  phase 

mycent roid_x=zeros ( [Nsub,  Nsub] ) ;  %initialise  centroid  values 
mycent roid_y=zeros ( [Nsub, Nsub] ) ; 

for  x_pos  =  l:Nsub  %perform  centroiding  for  each  subaperture 
for  y_pos  =  l:Nsub 

if  (my Tele scope .mask (x_pos, y_pos)>0) 
ellipse=myEllipseGen2 (myLaser,  . .  . 

1 ,  my Telescope .x (x_pos, y_pos) ,  .  . . 
myTele scope .y(x_pos,y_pos) ,  eta,  nu, . . . 
myLayer , minor_ax ,  1 ,  .  .  . 
screen_grad_x  (x_pos ,  y_pos )  ,  .  .  . 
screen_grad_y  (x_pos  ,  y_pos )  ,  shift fac)  ; 
%generates  beacon  image,  shifted  depending  on 
%turbulence  gradient 
convolved_ellipse=conv2 (PSF, ellipse) ; 
cutof f=round (0 . 25*length (convolved_ellipse) ) ; 
conv_ellipse4=convolved_ellipse  (cutoff :  .  .  . 

3*cutof f , cutoff : 3*cutof f ) ; 
integrated_convolved=myInt  ( conv_ellipse4 ,  .  .  . 
myTelescope . res ) ; 

sum_tot=sum  (sum  (integrated.convolved)  )  ; 
mycent roid_x  (x_pos  ,  y_pos  )  =sum  (sum  (  .  .  . 
integrated.convolved.  *  .  .  . 
myTelescope  .  x_loc/ sum.tot )  )  .  .  . 

/myTelescope . Diameter; 
mycent roid_y  (x_pos  ,  y_pos  )  =sum  (sum  (  .  .  . 
integrated_convolved .  *  .  .  . 
myTelescope  .  y_loc/sum_tot )  )  .  .  . 
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/myTelescope . Diameter; 


221 

222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 


else 

mycentroicLx  (x_pos  ,  y_pos  )  =0; 
mycentroid_y  (x_pos  ,  y_pos  )  =0; 

end 

end 

end 

mycent  roid_x_n=mycent  roid_x  .  / nocal_x2 ;  %apply  standard 
mycent  roid_y_n=mycentroid_y  .  / nocal_y2 ;  %calibration  factor 
mycent roid_x=mycentroid_x ./ cal_x; %apply  enhanced  calibration 
mycent roid_y=mycentroid_y .  /cal_y ;  %factor 

S_m=  [mycent roid_y  (ActiveSubs)  ;  mycentroid_x  (ActiveSubs)  ]  ; 
P_m=zeros (size (ActiveActs) )  ; 

P_m (mask) =H* S_m;  %reconstruct  enhaced  cal  phase 

S_n=  [mycent roid_y_n  (ActiveSubs)  ;  mycent roid_x_n  (ActiveSubs)  ]  ; 
P_n=zeros (size (ActiveActs) ) ; 

P_n (mask) =H* S_n ;  %reconstruct  standard  cal  phase 

tempi (: ,  :) =P— P_m;  %take  difference  between  real  and 

%enhanced  cal  phase  values 

temp2 ( : , : ) =P— P_n ;  %take  difference  between  real  and 

%standard  cal  phase  values 

EnhCal (idxprime, : ,  : ) =templ ;  %save  result  for  each  iteration 
StdCal (idxprime, : ,  :  ) =temp2; 

%results  still  need  analysis,  and  to  be  RMSE ' d 

end 

end 
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Listing  A. 2.  myEllipseGen.m 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 


function  [ellipse]  =  myEllipseGen (  Laser,  time,  x,  y, . . . 

eta,nu,  myLayer, minor.ax, A  ,  x.grad, y_grad, shift ) 
%MYELLIPSEGEN  generates  quassian  beacon  for  a  given  system 
%orientation . 

%  Eliiptical  shape  depends  on  the  location  of  the  subaperture 
%  and  laser,  local  tilt  on  the  wavefront,  speed  of  the  lenslet, 
%  sodium  layer  dimensions  and  look  angle 

o_ 

o 

%  @author:  FLTLT  Russell  McGuigan 
%  @date:  21Feb2011 

%  @ inputs: 

%  Laser:  Structure  containing  laser  Coude  rotation  info 

%  time:  position  along  track  (if  multiple  az/el  angles 

%  defined) ,  should  be  1  for  single  position 

%  x/y:  subaperture  geometry  [m] 

%  eta/nu:  grid  on  which  image  is  formed 

%  myLayer:  sodium  boundary  information 

%  minor.ax:  Ellipse  minor  axis  (fixed) 

%  A:  Ellipse  max  intensity  (usually  1) 

%  x_grad:  local  x  tilt  over  subaperture 

%  y_grad:  local  y  tilt  over  subaperture 

%  shift:  amount  that  tilt  affects  spot  displacement, 

%  similar  to  speed  or  focal  length  of  lenslet 


@outputs : 
ellipse : 


Gaussian  (elliptical  looking)  beacon  image 
which  depends  on  all  the  above  inputs 


theta_rel2  =  atan2  (Laser  .  beacon_y  (time) —y,  ..  . 

Laser . beacon_x (time) — x) ; 

%Calculate  distance  and  angle  to  laser  from  each  subaperture 
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32  dist_rel=sqrt ( (Laser . beacon_x (time)  —  x)  ,"2+.  .  . 

33  (Laser  .  beacon.y  (time)  — y)  .  ~2  )  ; 

34  %Calculate  ground  distance  from  subaperture  to 

35  %top  and  bottom  of  LGS 

36  theta_width=atan2  (myLayer  .  High  (time )  ,  dist_rel)  —  .  .  . 

37  atan2 (myLayer . Low (time) , dist.rel )  ; 

38  %calculate  angular  width  of  LGS  and  convert  to  distance 

39  theta2=theta_width*myLayer .Mean (time)  ; 

40  theta2=max  (theta2  ,  minor_ax)  ; 

41  %ensures  that  the  width  is  never  less  than  the  minor  axis 

42  %width 

43  x_of f set=x_grad*shift*cos  (theta_rel2)  — .  .  . 

44  y_grad*shift*sin (theta_rel2 ) ; 

45  %coordinate  transformation  to  rotate  ellipse  offset 

46  %to  align  with  beacon 

47  y_of f set=x_grad*shift*sin (theta_rel2 ) + . . . 

48  y_grad*shift*cos (theta_rel2 )  ; 

49  eta2=eta*cos (theta_rel2 ) — nu*sin (theta_rel2)  ; 
so  %coordinate  transformation  to  rotate  ellipse 

51  %to  align  with  beacon 

52  nu2=eta*sin (theta_rel2) +nu*cos (theta_rel2 ) ; 

53  ellipse  =  l*exp (— (  (eta2— x_of f set )  .  "2/ (2*theta2)  +  . . . 

54  (nu2— y_of  f  set )  .  "2/  (2*minor_ax)  )  )  *A; 

55  %generate  actual  ellipse 

56  end 


Listing  A. 3.  Miscellaneous  Functions 

1  function  [  myPath  ]  =  myTrack (  az,el  ) 

2  %MYTRACK  converts  azimuth  and  elevation  angle  data  into  a 

3  %single  structure  and  converts  from  degrees  to  radians 
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5  % 


@author : 


FLTLT  Russell  McGuigan 

6  %  @date:  21Feb2011 

7  %  @ inputs: 

8  %  az:  telescope  azimuthal  angle  [deg] 

9  %  el:  telescope  elevation  angle  [deg] 

10  %  @ outputs: 

11  %  myPath:  structure  storing  az  and  el  data  [rad] 

ir,  9'9'9'9'2'9'9'9'9'9'9'9'9'2'9'9'9'9'9'9'9'9'9'9'2'2'2'9'9'9'9'9'9'9'9'9'9'9'9'9'9' 

JLz  ooooooooooooooooooooooooooooooooooooooooo 

13  if  length (az ) ~=length (el ) ;  %basic  length  check  to  ensure 

14  %a  value  wasn't  missed 

15  disp (' Azimuth  and  Elevation  data  must  be  same  length') 

16  else 

17  for  idx=l :  length  (az  )  ; 

is  myPath . el ( idx) =el (idx) *pi/180 ; %convert  to  radians 

19  myPath . az ( idx) =az (idx) *pi/180 ; 

20  end 

21  end 

22 

23  end 

24 

25  function  [  Layers  ]  =  myLayers (  Path,  Low,  High,  wvl,  R  ) 

26  %MYLAYERS  Determines  high  low  and  mean  boundary  heights  for 

27  %a  given  satellite  look  angle 

28  % 


29 

o, 

o 

Sauthor : 

FLTLT  Russell 

McGuigan 

30 

o. 

o 

(Mate : 

21Feb2011 

31 

o_ 

o 

@ inputs : 

32 

o, 

o 

Path : 

telescope  track  data  (az  and 

el)  [rad] 

33 

o, 

o 

Low : 

sodium  layer 

lower  boundary 

at  zenith 

[m] 

34 

o_ 

o 

High: 

sodium  layer 

upper  boundary 

at  zenith 

[m] 

35 

o. 

o 

wvl : 

sodium  laser 

wavelength  [m] 
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36  % 


R: 


Earth  radius  [m] 


37  %  @ outputs: 

38  %  Layers:  structure  containing  boundary  height  info 

o9  ooooooooooooooooooooooooooooooooooooooooo 

40 

41  if  nargin  <5  %generic  earth  radius  if  not  specified 

42  R=  6371000; 

43  end 

44  if  nargin<4  %sodium  wavelength  (if  not  specified) 

45  wvl=58 9 . 2e— 9; 

46  end 

47  if  nargin==l;  %zenith  boundary  heights  (if  not  specified) 

48  Low=85000; 

49  High=95000; 
so  end 

si  Layers. R=R;  %stores  earth  radius  and  laser  wavelength 

52  Layers  .  wvl=wvl ; 

53  Layers . High=Layers .R*cos(pi/2+Path.el)+. . . 

54  sqrt ( (High+Layers . R) "2— (Layers . R*sin (pi/2+Path . el ) )  . ~  2 ) ; 

55  %Distance  to  top  boundary  layer,  assuming 

56  %spherical  earth 

57  Layers . Low=Layers . R*cos (pi/ 2+Path . el ) + . . . 

58  sqrt ( (Low+Layers . R) ~2— (Layers . R*sin(pi/2+Path.el) ) ."2) ; 

59  %Distance  to  lower  boundary  layer 

60  Layers .Mean=0 . 5* (Layers . High+Layers . Low) ; %Layer  midpoint 

61  end 

62 

63  function  [Telescope] =myTel (Nsub, Diameter, L_sep , FOV, res, theta_0 ) 

64  %MYTEL  Stores  telescope  data  into  a  single  structure  and 

65  %generates  detector  geometry 

66  % 


67  %  @author:  FLTLT  Russell  McGuigan 


68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 


@date:  21Feb2011 

@inputs : 

Nsub:  Number  of  telescope  subapertures 

Diameter:  Telescope  diameter  [m] 

L_sep:  radial  separation  of  laser  from  edge  of 

telescope  [m] .  Put  —Diameter/2  for  centre 
launched 

FOV :  Telescope  FOV  (no  longer  used) 

res:  detector  resolution.  2  for  quad  cell 

theta_0 :  initial  Coude  offset  of  exit  pupil  [deg] 

@ outputs : 

Telescope:  structure  containing  telescope  information 

ans  subaperture  and  detector  geometries 


Telescope .Nsub=Nsub;  %stores  data  into  structure 
Telescope . Diameter=Diameter ; 

Telescope .  L_sep=L_sep; 

Telescope . F0V=F0V; 

Telescope . res=res ; 

Telescope . theta_0=theta_0 *pi/180 ; %convert  offset  to  radians 
x= (—Nsub/ 2+1/ 2 : 1 : Nsub/ 2— 1/2) *Diameter/Nsub; 

%generate  subaperture  geometry 
[x  y ] =meshgrid (x) ; 

Telescope . x=x; 

Telescope . y=y; 

Telescope .mask=double (x . ~2+y . ~2<= (Diameter/ 2 ) "2 ) ; 

%generate  subaperture  mask 

x_loc=Diameter/2* (  (—  (res— 1) /res )  : 2/res :  ( (res— 1) /res )  )  ; 

%generate  detector  geometry 
[Telescope . x_loc  Telescope . y_loc]  =  meshgrid (x_loc) ; 
end 
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100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 

125 

126 

127 

128 

129 

130 

131 


function  [  atmosphere  ]  =  my Atmosphere (Nsub, Diameter , rO ,  N) 
%MYATMOSPHERE  generates  random  phase  screen,  calculates  its 
%gradient . 

o. 

o 

%  @author:  FLTLT  Russell  McGuigan 
%  @date:  21Feb2011 

%  @ inputs: 

Number  of  telescope  subapertures 
Telescope  diameter  [m] 

Atmospheric  coherence  diamter  (Fried's 
parameter  [m] 

N:  Number  of  grid  points  across  phase  screen 

@outputs : 

atmosphere:  structure  containing  phase  screen  values 

and  average  gradients  over  each  subaperture 


Nsub : 
Diameter : 
rO  : 


PHZ=kolmogphz screen (N, Diameter,  rO ,  0 ,  [ 0  0 ] , randi ( 10000 ) ) ; 

%generates  phase  screen  for  given  parameters,  no  added 
%tilt 

phz=PHZ .phase; 

steps=round ( 1 inspace ( 1 , length (phz ) , Nsub+1 ) ) ; 

%determines  locations  over  which  to  average  gradients 
atmosphere  .  screen_grad_x=zeros  (Nsub)  ;  %initialis  gradients 
atmosphere  .  screen_grad_y=atmosphere  .  screen_grad_x; 
for  idxl=l:Nsub  %loop  for  each  subaperture 
for  idx2=l:Nsub 

sub_screen=phz (steps (idxl) : steps (idxl+1) , . . . 
steps (idx2) : steps (idx2+l) ) ; 

%sub  divide  phase  screen  into  subapertures 
[ fx,  fy ] =gradient ( sub-screen ) ;  %compute  gradients  for 
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132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 


%each  subaperture 

atmosphere  .  screen_grad_x  ( idxl ,  idx2 )  =mean  (mean  ( f x)  )  ; 

%average  the  gradient  over  the  subaperture  (x  and  y) 
atmosphere  .  s creen.gr ad_y  ( idxl ,  idx2 )  =mean  (mean  ( fy )  )  ; 
end 

end 

atmosphere . phase=phz ;  %store  phase  values 

end 


function  [  PSF  eta  nu]  =  myPSFgen (  Telescope,  res, a  ) 
%MYPSFGEN  generates  diffraction  limited  spot 

o, 

o 

%  Sauthor:  FLTLT  Russell  McGuigan 
%  @date:  21Feb2011 

%  @ inputs: 

%  Telescope:  Structure  containing  telescope  data 
%  res:  number  of  grid  points  for  PSF  generation 

%  a:  width  of  PSF 

%  @outputs: 

%  PSF:  sinc~2  diffraction  limited 

%  eta:  cartesian  grid  of  x  values  on  which  the 

%  PSF  lies 

%  nu :  eta's  y  value  equivalent 


eta=linspace (— 1 . 75, 1 . 75, res) ;  %initialise  sensor  grid 
[eta  nu]  =  meshgrid (eta) ; 

PSF=sinc (a*eta) ,"2.*sinc(a*nu) ."2; 

%generate  PSF  (sinc~2)  function 

end 

function  [  Laser  ]  =  myCLrotate (  Telescope,  Track) 
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164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 

183 


%MYC_ROTATE  determines  the  Coude  rotation  of  the  laser 


%beacon  and  it 1 s  x  and  y  position 

o_ 

o 

%  @author:  FLTLT  Russell  McGuigan 

%  Sdate:  21Feb2011 

%  @ inputs: 

%  Telescope:  Structure  containing  telescope  data 

%  Track:  Structure  containing  target  azimuth 

%  and  elvation  values  in  radians 

%  Soutputs: 

%  Laser:  Coude  angle  and  position  of  laser 

9'2-2'9'2'9-9'9'9'9'2-2'2-2'9'9'9-9'9'9'2'2-9'9'2'9'9'9'9'2'9'9'2-9-2'9-9'9'9'9'2' 

ooooooooooooooooooooooooooooooooooooooooo 


Laser . theta_beac=Track . a  z— Track . el+Tele scope . theta_0 ; 
%determine  Coude  rotation 

Laser . beacon_x= (Telescope . Diameter/2  +  Telescope  .  L_sep)  . . . 

*cos (Laser . theta_beac) ; %determine  laser  x  location 
Laser . beacon_y= (Telescope . Diameter/2  +  Telescope  .  L_sep)  . . . 
*sin (Laser . theta_beac)  ;  %determine  laser  y  location 

end 
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